This is a preprint of: “Global phase portraits of the generalized van der Pol systems”, Jaume Llibre, Clàudia Valls, Bull. Sci. Math., vol. 182, 103213, 2023. DOI: [10.1016/j.bulsci.2022.103213] GLOBAL PHASE PORTRAITS OF THE GENERALIZED VAN DER POL SYSTEMS JAUME LLIBRE1 AND CLAUDIA VALLS2 Abstract. We consider the generalized van der Pol systems ẋ = y, ẏ = −x + (1 − x2 )f (y), where f ∈ R[y]. The classical van der Pol systems have f (y) = y. We first characterize when the origin of the generalized van der Pol systems is a center, and second we provide the global phase portraits in the Poincaré disc of the generalized van der Pol when f (y) = a1 y + a2 y 2 for all a1 , a2 ∈ R. 1. Introduction and statement of the main results In this paper we deal with the generalized van der Pol systems (1) ẋ = y, ẏ = −x + (1 − x2 )f (y), Pn i where f (y) is the polynomial i=1 ai y with n ≥ 1. The classical van der Pol equation has f (y) = a1 y. There is no doubt about the importance of this differential system and this is one of the reasons why it has been studied for so many authors. For instance, if one enters the four words van der Pol, differential, equation or system in MathSciNet, one would receive 768 articles at the time that paper is being written. For instance some four recent papers on variations on the van der Pol system are [3, 7, 8, 9]. The main two theorems of this paper are the following ones. Theorem 1. System (1) has a center at the origin if and only if the polynomial f (y) is even. The proof of Theorem 1 will be given in section 3. Theorem 2. The global phase portrait of system (1) with f (y) = a1 y + a2 y 2 with a1 , a2 ∈ R is topologically equivalent to the one of : (i) Figure 1(a) if a2 = 0 and a1 ̸= 0; (ii) Figure 1(b) if a2 ̸= 0 and a1 = 0. Moreover, there are systems (1) with a2 a1 ̸= 0 whose global phase portrait is topo- logically equivalent to the one of Figure 1(c). 2010 Mathematics Subject Classification. 34C05. Key words and phrases. van der Pol systems, center, Lyapunov constant, Poincaré compacti- fication, blow ups. 1 2 J. LLIBRE AND C. VALLS (a) a2 = 0, a1 ̸= 0 (b) a2 ̸= 0, a1 = 0 (c) a2 a1 ̸= 0 Figure 1. The three phase portraits in the Poincaré disc of Theorem 2. The proof of Theorem 2 will be given in section 4. It was proved in [6, 10] that system (1) with a2 a1 = 0 has always a unique limit cycle. Based in numerical evidence it seems that system (1) with a2 a1 ̸= 0 has also a unique limit cycle, so we conjecture that system (1) has always a unique limit cycle which allows us to state the following conjecture on the global phase portraits of system (1) with a2 a1 ̸= 0. Conjecture 3. The global phase portrait of system (1) with f (y) = a1 y + a2 y 2 and a2 a1 ̸= 0 is topologically equivalent to the one of Figure 1(c). Before proving Theorems 1 and 2 we have a preliminary section, section 2, where we have the notions and results for proving such theorems. 2. Preliminary results 2.1. Singular points. Consider a differential system of the form (2) ẋ = p(x, y), ẏ = q(x, y) being p, q ∈ R[x, y]. The point (a, b) is a singular point of the differential system (2) if p(a, b) = q(a, b) = 0. The singular point (a, b) is hyperbolic if the eigenvalues of the Jacobian matrix of the function (p, q) evaluated at (a, b) have non-zero real part. The classification of the local phase portraits of the hyperbolic singular points is well known, see for instance [4, Theorem 2.15]. In this paper when we characterize the local phase portrait of a hyperbolic singular point we will use that theorem. The singular point (a, b) is semi-hyperbolic if one and only one of the eigenvalues of the Jacobian matrix of the function (p, q) evaluated at (a, b) is zero. Also the classification of the local phase portraits of the semi-hyperbolic singular points is well known, see for instance [4, Theorem 2.19]. Again in this paper when we characterize the local phase portrait of a semi-hyperbolic singular point we will use that theorem. Consider the differential system X∞ ∞ X (3) ẋ = pi (x, y), ẏ = qi (x, y), i=1 i=1 GLOBAL PHASE PORTRAITS OF THE GENERALIZED VAN DER POL SYSTEMS 3 where pi and qi are homogeneous polynomials of degree i, for i ≥ 1. The character- istic directions of the singular point localized at the origin of coordinates of system (3) are given by the straight lines trough the origin defined by the real linear factors of the homogeneous polynomial pk (x, y)y − qk (x, y)x, where k is the minimum i for which the polynomials pi or qi are non-zero. It is known that the orbits which end or start at the origin of coordinates must arrive or exit tangent to these straight lines. For more details on the characteristic directions see for example [2]. When the Jacobian matrix of the function (p, q) evaluated at the singular point (a, b) of the differential system (2) is identically zero, then the singular point is called linearly zero and the local phase portrait at this singular point can be studied doing special changes of variables called blow ups, see for instance [1]. Here we use vertical blow ups and when the vertical direction is a characteristic direction we twist it to another direction in order that in the new variables the vertical axis is not a characteristic direction. Let Φt be a smooth flow on a manifold M and let C be a submanifold of M consisting entirely of singular points of the flow. C is called normally hyperbolic if the tangent bundle to M over C splits into three subbundles T C, E s and E u invariant under the differential dΦt and satisfying (i) dΦt contracts E s exponentially, (ii) dΦt expandsE u exponentially, (iii) T C is the tangent bundle of C. For normally hyperbolic submanifolds one has the usual existence of smooth stable and unstable manifolds together with the persistence of these invariant manifolds under small perturbations. More precisely, we have the following theorem, for a proof see [5]. Theorem 4. Let C be a normally hyperbolic submanifold of singular points for the flow Φt . Then there exist smooth stable and unstable manifolds tangent along C to E s ⊕ T C and E u ⊕ T C, respectively. Moreover, both C and the stable and unstable manifolds are persistent under small perturbations of the flow. 2.2. The Poincaré compactification. Roughly speaking the Poincaré compact- ification consists in identifying the plane R2 with the interior of a closed unit disc centered at the origin of coordinates, called the Poincaré disc. Then the boundary of this disc (the unit circle centered at the origin) is identified with the infinity of R2 . Note that in R2 we can go or come from the infinity in as many as directions as points has that circle. In order to classify the global dynamics of a polynomial differential system one of the main steps is to characterize the local phase portraits of its finite and infinite singular points in the Poincaré disc. For doing this we need the equations of our polynomial differential systems initially in R2 in the Poincaré disc. Consider the differential system (2) in R2 , where p and q are real polynomials in the variables x and y of degrees d1 and d2 , respectively. Then the degree of the polynomial differential system (2) is d = max{d1 , d2 }. 4 J. LLIBRE AND C. VALLS Denote by Tp S2 be the tangent space to the 2-dimensional sphere S2 = {s = (s1 , s2 , s3 ) ∈ R3 : s21 + s22 + s23 = 1} at the point p, we call this sphere the Poincaré sphere. We consider that the polynomial differential system (2) is defined in the tangent plane to S2 at the point (0, 0, 1), i.e. we have identified R2 with T(0,0,1) S2 . The central projection f : T(0,0,1) S2 → S2 send each point p of T(0,0,1) S2 to two points of S2 , one in the northern hemisphere and the other in the southern hemisphere. These two points are the intersection of the straight line through p and the origin of coordinates (the center of the sphere). So the map f defines two copies of the polynomial differential system (2) on the sphere, one in the open northern hemisphere and the other in the open southern hemisphere. If X = (p, q) is the vector field associated to the polynomial differential system (2), we denote by X ′ the vector field Df ◦ X defined on S2 except on its equator S1 = {s ∈ S2 : s3 = 0}. Clearly S1 can be identified with the infinity of R2 . If the degree of the polynomial vector fiedl X is d, then p(X) is the only analytic extension of sd−1 3 X ′ to S2 . The vector field p(X) on S2 is called the Poincaré compactification of the vector field X, for more details see [4, chapter 5]. On the Poincaré sphere S2 we use the following six local charts, which are given by Ui = {s ∈ S2 : si > 0} and Vi = {s ∈ S2 : si < 0}, for i = 1, 2, 3, with the corresponding diffeomorphisms φi : Ui → R2 , ψi : Vi → R2 , defined by φi (s) = ψi (s) = (sm /si , sn /si ) = (u, v) for m < n and m, n ̸= i. Thus the coordinates (u, v) will play different roles in the distinct local charts. The expressions of the vector field p(X) are         1 u 1 u 1 u (u̇, v̇) = v d Q , − uP , , −v d+1 P , in U1 , v v v v v v         u 1 u 1 u 1 (u̇, v̇) = v d P , − uQ , , −v d+1 Q , in U2 , v v v v v v (u̇, v̇) = (P (u, v), Q(u, v)) in U3 . We note that the expressions of the vector field p(X) in the local chart (Vi , ψi ) is equal to the expression in the local chart (Ui , ϕi ) multiplied by (−1)d−1 for i = 1, 2, 3. The orthogonal projection under π(y1 , y2 , y3 ) = (y1 , y2 ) of the closed northern hemisphere of S2 onto the plane s3 = 0 is a closed disc D2 of radius one centered at the origin of coordinates called the Poincaré disc. Since a copy of the vector field X on the plane R2 is in the open northern hemisphere of S2 , the interior of the Poincaré disc D2 is identified with R2 and the boundary of D2 , the equator S1 of S2 , is identified with the infinity of R2 . Consequently the phase portrait of the vector field X extended to the infinity corresponds to the projection of the phase portrait of the vector field p(X) on the Poincaré disc D2 . The singular points of p(X) in the Poincaré disc lying on S1 are the infinite singular points of the vector field X. The singular points of p(X) in the interior of the Poincaré disc, i.e. on D2 \ S1 , are the finite singular points. We note that in GLOBAL PHASE PORTRAITS OF THE GENERALIZED VAN DER POL SYSTEMS 5 the local charts U1 , U2 , V1 and V2 the infinite singular points have their coordinate v = 0. For a polynomial differential system (2) if s ∈ S1 is an infinite singular point, then −s ∈ S1 is another infinite singular point. Thus the number of infinite singular points is even. 2.3. The Poincaré-Bendixson Theorem. In what follows we are going to as- sume that ∆ is an open subset of R2 and X is a vector field of class C r with r ≥ 1. Also, in ∆, γp+ denotes a positive semi–orbit passing through the point p. Let φ(t) = φ(t, p) = φp (t) be the integral curve of X passing through the point p, defined on its maximal interval Ip = (α, ω). If ω = ∞ we define the set ω(p) = {q ∈ ∆ : there exist {tn } with tn → ∞ and φ(tn ) → q when n → ∞}. In the same way, if α = −∞ we define the set α(p) = {q ∈ ∆ : there exist {tn } with tn → −∞ and φ(tn ) → q when n → ∞}. The sets ω(p) and α(p) are called the ω–limit set and the α–limit set of p, respec- tively. Theorem 5 (Poincaré–Bendixson Theorem). Let φ(t) = φ(t, p) be an orbit of X defined for all t ≥ 0, such that γp+ is contained in a compact set K ⊂ ∆. Assume that the vector field X has at most a finite number of singularities in K. Then one of the following statements holds. (i) If ω(p) contains only regular points, then ω(p) is a periodic orbit. (ii) If ω(p) contains both regular and singular points, then ω(p) is formed by a set of orbits, every one of which tends to one of the singular points in ω(p) as t → ±∞. (iii) If ω(p) does not contain regular points, then ω(p) is a unique singular point. For a proof of the Poincaré–Bendixson Theorem see for instance [4, Theorem 1.25]. 2.4. Uniqueness of limit cycles. The next result is proved in [6], see also [10, Theorem 4.1]. Theorem 6. Consider the differential system (4) ẋ = y, ẏ = −g(x) − f (x)y. Z x Let F (x) = f (s)ds. Assume that the following conditions hold. 0 (i) f : R → R is continuous and there exist a < 0 < b such that f (x) < 0 if x ∈ (a, b), and f (x) > 0 if x < a and x > b; (ii) there exists c > 0 such that F (c) = F (−c) = 0; (iii) F (+∞) = +∞, or F (−∞) = −∞. Then the differential system (4) has a unique limit cycles which is stable. 6 J. LLIBRE AND C. VALLS 3. Proof of Theorem 1 Assume that the polynomial f (y) is even. Then the differential system (1) is invariant under the change of variables (x, y, t) 7→ (x, −y, −t). Therefore, since the origin is monodromic in a neighborhood of it, the origin must be a center. Now assume that the polynomial f (y) is not even, so we can write it as f (y) = fe (y)+fo (y), where fe (y) and fo (y) are the even and the odd part of the polynomial f (y), respectively. Consider the determinant y −x + fe (y)(1 − x2 ) = yfo (y)(1 − x2 ). y −x + (fe (y) + fo (y))(1 − x2 ) Since yfo (y)(1 − x2 ) has a constant sign in a convenient sufficiently small neighbor- hood U of the origen, at any point of U the oriented angle between the vector field (y, −x + fe (y)(1 − x2 )) having a center in U and the vector field (y, −x + (fe (y) + fo (y))(1 − x2 )) has the same sign. Consequently the vector field (y, −x + (fe (y) + fo (y))(1 − x2 )) has an unstable focus at the origin if yfo (y)(1 − x2 ) < 0 in U , and it has a stable focus at the origin if yfo (y)(1 − x2 ) > 0. Consequently the vector field has a center at the origin if and only if fo (y) = 0. This completes the proof of Theorem 1. 4. Proof of Theorem 2 The unique finite singular point of system (1) is the origin. If a1 ∈ (−2, 2) it is a strong focus if a1 ̸= 0, and if a1 = 0 in view of Theorem 1 it is a center. Moreover if |a2 | > 2 then the origin is a node. More precisely, the origin is a stable node if a1 ∈ (−∞, 2], a stable focus if a1 ∈ (−2, 0), a center if a1 = 0, an unstable focus if a1 ∈ (0, 2) and an unstable node if a1 ∈ [2, +∞). Now we study the infinite singular points for each of the cases in the statements of the theorem. 4.1. The infinite singular points of the local charts U1 and U2 when a2 = 0 and a1 ̸= 0. Now we shall study the infinite singular points of the differential system (1) in the local charts U1 and U2 . We note that in this case doing the change of variables (x, y, t) → (x, −y, −t) we can assume without loss of generality that a1 > 0. 4.1.1. The local chart U1 . System (1) in the chart U1 writes (5) u̇ = −a1 u − v 2 + a1 uv 2 − u2 v 2 , v̇ = −uv 3 . The unique infinite singular point in this chart is the origin. The eigenvalues of the linear part of system (5) at the singular point (0, 0) are 0 and −a1 and so it is semi-hyperbolic singular point. Applying [4, Theorem 2.19] we obtain that it is a semi-hyperbolic saddle. GLOBAL PHASE PORTRAITS OF THE GENERALIZED VAN DER POL SYSTEMS 7 4.1.2. The local chart U2 . System (1) in the chart U2 writes (6) u̇ = v 2 + a1 u3 − a1 uv 2 + u2 v 2 , v̇ = v(a1 u2 − a1 v 2 + uv 2 ). The origin is a singular point. The linear part of system (6) at the origin is iden- tically zero, so in order to determine its local phase portrait we must do blow ups. The characteristic directions at the origin of U2 are given by the real linear factors of v 3 . Since u = 0 is not a characteristic direction we do the vertical blow up (u, v) → (u1 , u1 v1 ) and system (6) writes in the new variables u̇1 = −u21 (−a1 u1 − v12 + a1 u1 v12 − u21 v12 ), v̇1 = −u1 v13 . Doing a rescaling of the independent variable we eliminate the common factor u1 between u̇1 and v̇1 and we obtain the system u̇1 = −u1 (−a1 u1 − v12 + a1 u1 v12 − u21 v12 ), (7) v̇1 = −v13 . The unique singular point of system (7) on the straight line u1 = 0 is (0, 0). The linear part of system (7) evaluated at (0, 0) is linearly zero so we must do another blow up. Then the vertical axis u1 = 0 is a characteristic direction at the origin of the local chart U2 . Therefore before doing a vertical blow up we translate the direction u1 = 0 to u1 = v1 doing the change of variables (u1 , v1 ) = (u2 − v2 , v2 ). Then system (7) becomes (8) u̇2 = a1 u22 − 2a1 u2 v2 + a1 v22 + u2 v22 − a1 u22 v22 + u32 v22 − 2v23 + 2a1 u2 v23 − 3u22 v23 −a1 v24 + 3u2 v24 − v25 , v̇2 = −v23 . Now we do the vertical blow up (u2 , v2 ) → (u3 , u3 v3 ) and system (8) writes in the new variables u̇3 = −u23 (−a1 + 2a1 v3 − a1 v32 − u3 v32 + a1 u23 v32 − u33 v32 + 2u3 v33 − 2a1 u23 v33 +3u33 v33 + a1 u23 v34 − 3u33 v34 + u33 v35 ), v̇3 = u3 (v3 − 1)v3 (a1 − a1 v3 + 2u3 v32 − a1 u23 v32 + u33 v32 + a1 u23 v33 − 2u33 v33 + u33 v34 ). Doing a rescaling of the time we eliminate the common factor u3 between (u̇3 , v̇3 ) and we get the system (9) u̇3 = −u3 (−a1 + 2a1 v3 − a1 v32 − u3 v32 + a1 u23 v32 − u33 v32 + 2u3 v33 − 2a1 u23 v33 +3u33 v33 + a1 u23 v34 − 3u33 v34 + u33 v35 ), v̇3 = (v3 − 1)v3 (a1 − a1 v3 + 2u3 v32 − a1 u23 v32 + u33 v32 + a1 u23 v33 − 2u33 v33 +u33 v34 ). The singular points of system (9) on the straight line u3 = 0 are (0, 0) and (0, 1). The eigenvalues of the linear part of system (9) evaluated at the singular point (0, 0) are −a1 and a1 , and so it is a saddle. The linear part of system (9) evaluated at (0, 1) is linearly zero so we must do another blow-up. First we translate the singular point (0, 1) to the origin of coordinates in order to study its local phase portrait 8 J. LLIBRE AND C. VALLS doing the change (u3 , v3 ) = (u4 , 1 + v4 ). Therefore system (9) in the variables (u4 , v4 ) becomes (10) u̇4 = −u4 (u4 + 4u4 v4 − a1 v42 + 5u4 v42 + a1 u24 v42 + 2u4 v43 + 2a1 u24 v43 + u34 v43 +a1 u24 v44 + 2u34 v44 + u34 v45 ), v̇4 = 2u4 v4 − a1 v42 + 6u4 v42 + a1 u24 v42 − a1 v43 + 6u4 v43 + 3a1 u24 v43 + u34 v43 + 2u4 v44 +3a1 u24 v44 + 3u34 v44 + a1 u24 v45 + 3u34 v45 + u34 v46 . The characteristic directions at the origin are the real linear factors of u4 v4 (a1 v4 − 3u4 ) = 0. Then the vertical axis u4 = 0 is a characteristic direction at the origin of system (10). Before doing a vertical blow up we translate the direction u4 = 0 to u4 = v4 doing the change of variables (u4 , v4 ) = (u5 − v5 , v5 ). Then system (10) becomes (11) u̇5 = −u25 + 4u5 v5 − 4u25 v5 − 3v52 − a1 v52 + 14u5 v52 + a1 u5 v52 − 5u25 v52 + a1 u25 v52 −a1 u35 v52 − 10v53 − 2a1 v53 + 16u5 v53 − 2a1 u5 v53 − 2u25 v53 + 6a1 u25 v53 + u35 v53 −2a1 u35 v53 − u45 v53 − 11v54 + a1 v54 + 6u5 v54 − 9a1 u5 v54 − 3u25 v54 + 9a1 u25 v54 +7u35 v54 − a1 u35 v54 − 2u45 v54 − 4v55 + 4a1 v55 + 3u5 v55 − 12a1 u5 v55 − 15u25 v55 +4a1 u25 v55 + 11u35 v55 − u45 v55 − v56 + 5a1 v56 + 13u5 v56 − 5a1 u5 v56 − 21u25 v56 +5u35 v56 − 4v57 + 2a1 v57 + 17u5 v57 − 9u25 v57 − 5v58 + 7u5 v58 − 2v59 , v̇5 = 2u5 v5 − 2v52 − a1 v52 + 6u5 v52 + a1 u25 v52 − 6v53 − a1 v53 + 6u5 v53 − 2a1 u5 v53 +3a1 u25 v53 + u35 v53 − 6v54 + a1 v54 + 2u5 v54 − 6a1 u5 v54 − 3u25 v54 + 3a1 u25 v54 +3u35 v54 − 2v55 + 3a1 v55 + 3u5 v55 − 6a1 u5 v55 − 9u25 v55 + a1 u25 v55 + 3u35 v55 − v56 +3a1 v56 + 9u5 v56 − 2a1 u5 v56 − 9u25 v56 + u35 v56 − 3v57 + a1 v57 + 9u5 v57 − 3u25 v57 −3v58 + 3u5 v58 − v59 . Now we do the vertical blow up (u5 , v5 ) → (u6 , u6 v6 ) and system (11) writes u̇6 = −u26 1 − 4v6 + 4u6 v6 + 3v62 + a1 v62 − 14u6 v62 − a1 u6 v62 + 5u26 v62 − a1 u26 v62 +a1 u36 v62 + 10u6 v63 + 2a1 u6 v63 − 16u26 v63 + 2a1 u26 v63 + 2u36 v63 − 6a1 u36 v63 − u46 v63 +2a1 u46 v63 + u56 v63 + 11u26 v64 − a1 u26 v64 − 6u36 v64 + 9a1 u36 v64 + 3u46 v64 − 9a1 u46 v64 −7u56 v64 + a1 u56 v64 + 2u66 v64 + 4u36 v65 − 4a1 u36 v65 − 3u46 v65 + 12a1 u46 v65 + 15u56 v65 −4a1 u56 v65 − 11u66 v65 + u76 v65 + u46 v66 − 5a1 u46 v66 − 13u56 v66 + 5a1 u56 v66 + 21u66 v66  −5u76 v66 + 4u56 v67 − 2a1 u56 v67 − 17u66 v67 + 9u76 v67 + 5u66 v68 − 7u76 v68 + 2u76 v69 , v̇6 = u6 (v6 − 1)v6 − 3 + 3v6 + a1 v6 − 10u6 v6 − a1 u26 v6 + 10u6 v62 + 2a1 u6 v62 −11u26 v62 + 2a1 u26 v62 − 4a1 u36 v62 − u46 v62 + 11u26 v63 − a1 u26 v63 − 4u36 v63 + 8a1 u36 v63 +3u46 v63 − 5a1 u46 v63 − 4u56 v63 + 4u36 v64 − 4a1 u36 v64 − 3u46 v64 + 10a1 u46 v64 + 12u56 v64 −2a1 u56 v64 − 5u66 v64 + u46 v65 − 5a1 u46 v65 − 12u56 v65 + 4a1 u56 v65 + 15u66 v65 − 2u76 v65  +4u56 v66 − 2a1 u56 v66 − 15u66 v66 + 6u76 v66 + 5u66 v67 − 6u76 v67 + 2u76 v68 . GLOBAL PHASE PORTRAITS OF THE GENERALIZED VAN DER POL SYSTEMS 9 Eliminating the common factor u6 between u̇6 and v̇6 rescaling the time we obtain the system (12) u̇6 = −u6 1 − 4v6 + 4u6 v6 + 3v62 + a1 v62 − 14u6 v62 − a1 u6 v62 + 5u26 v62 − a1 u26 v62 +a1 u36 v62 + 10u6 v63 + 2a1 u6 v63 − 16u26 v63 + 2a1 u26 v63 + 2u36 v63 − 6a1 u36 v63 − u46 v63 +2a1 u46 v63 + u56 v63 + 11u26 v64 − a1 u26 v64 − 6u36 v64 + 9a1 u36 v64 + 3u46 v64 − 9a1 u46 v64 −7u56 v64 + a1 u56 v64 + 2u66 v64 + 4u36 v65 − 4a1 u36 v65 − 3u46 v65 + 12a1 u46 v65 + 15u56 v65 −4a1 u56 v65 − 11u66 v65 + u76 v65 + u46 v66 − 5a1 u46 v66 − 13u56 v66 + 5a1 u56 v66 + 21u66 v66  −5u76 v66 + 4u56 v67 − 2a1 u56 v67 − 17u66 v67 + 9u76 v67 + 5u66 v68 − 7u76 v68 + 2u76 v69 , v̇6 = (v6 − 1)v6 − 3 + 3v6 + a1 v6 − 10u6 v6 − a1 u26 v6 + 10u6 v62 + 2a1 u6 v62 −11u26 v62 + 2a1 u26 v62 − 4a1 u36 v62 − u46 v62 + 11u26 v63 − a1 u26 v63 − 4u36 v63 + 8a1 u36 v63 +3u46 v63 − 5a1 u46 v63 − 4u56 v63 + 4u36 v64 − 4a1 u36 v64 − 3u46 v64 + 10a1 u46 v64 + 12u56 v64 −2a1 u56 v64 − 5u66 v64 + u46 v65 − 5a1 u46 v65 − 12u56 v65 + 4a1 u56 v65 + 15u66 v65 − 2u76 v65  +4u56 v66 − 2a1 u56 v66 − 15u66 v66 + 6u76 v66 + 5u66 v67 − 6u76 v67 + 2u76 v68 . The unique singular points of system (12) on the straight line u6 = 0 are (0, 0), (0, 1) and (0, 3/(a1 +3)). The eigenvalues of the linear part of system (12) evaluated at the singular point (0, 0) are 3 and −1, and so it is a hyperbolic saddle. Moreover the eigenvalues of the linear part of system (12) evaluated at the singular point (0, 1) are −a1 and a1 , and so it is also a hyperbolic saddle. Finally the eigenvalues of the linear part of system (12) evaluated at the singular point (0, 3/(a1 + 3)) are −3a1 /(3+a1 ) and −a1 /(3+a1 ) and so it is a stable node. Therefore the local phase portrait near the straight line u6 = 0 for system (12) is topologically equivalent to the one of Figure 2(a). Undoing the rescaling dt3 = u6 dt2 we get the phase portrait of Figure 2(b). Going back through the changes of variables from the phase portrait of Figure 2(b), we obtain the local phase portrait at the origin of system (11) which is topologically equivalent to the one of Figure 2(c). Going back through the changes of variables from the phase portrait of Fig- ure 2(c), we obtain the local phase portrait at the origin of system (10) which is topologically equivalent to the one of Figure 2(d). Going back through the changes of variables from the phase portrait of Figure 2(d), we obtain the local phase portrait around the straight line u3 = 0 which is topologically equivalent to the one of Figure 2(e). Again undoing the rescaling dt2 = u3 dt1 we obtain the phase portrait of Fig- ure 2(f). Going back through the changes of variables from the phase portrait of Figure 2(f), we obtain the local phase portrait at he origin of system (8) which is topologically equivalent to the one of Figure 2(g). Going back through the changes of variables from the phase portrait of Fig- ure 2(g), we obtain the local phase portrait at he origin of system (7) which is topologically equivalent to the one of Figure 2(h). 10 J. LLIBRE AND C. VALLS Undoing the rescaling dt1 = u1 dt we get the phase portrait of Figure 2(i). Finally going back through the changes of variables from the phase portrait of Figure 2(i), we obtain the local phase portrait at the origin of system (6) which is topologically equivalent to the one of Figure 2(j). Hence the origin of the local chart U2 has a local phase portrait which is topologically equivalent to an unstable node. (a) (b) (c) (d) (e) (f) (g) (h) (i) (j) Figure 2. Figures of the blow up of the singular point located at the origin of the local chart U2 of system (5). 4.2. The infinite singular points of the local charts U1 and U2 when a2 ̸= 0 and a1 = 0. We shall study the infinite singular points of the differential system (1) in the local charts U1 and U2 . Note that in this case doing the change of variables (x, y) → (−x, −y) we can assume, without loss of generality, that a2 > 0. 4.2.1. The local chart U1 . System (1) in the chart U1 writes (13) u̇ = −v 3 + u2 (−a2 + a2 v 2 − v 3 ), v̇ = −uv 4 GLOBAL PHASE PORTRAITS OF THE GENERALIZED VAN DER POL SYSTEMS 11 The unique infinite singular point in this chart is the origin. The linear part of system (13) at the origin is identically zero, so in order to determine its local phase portrait we must do blow ups. The characteristic directions at the origin of U2 are given by the real linear factors of u2 v = 0. Since u = 0 is a characteristic direction before doing a vertical blow up we translate the direction u = 0 to u = v doing the change of variables (u, v) = (u1 − v1 , v1 ). Then system (13) becomes (14) u̇1 = −a2 u21 + 2a2 u1 v1 − a2 v12 + a2 u21 v12 − v13 − 2a2 u1 v13 − u21 v13 + a2 v14 + u1 v14 , v̇1 = −u1 v14 + v15 . Now we do the vertical blow up (u1 , v1 ) → (u2 , u2 v2 ) and system (14) writes in the new variables u̇2 = u22 (−a2 + 2a2 v2 − a2 v22 + a2 u22 v22 − u2 v23 − 2a2 u22 v23 − u32 v23 + a2 u22 v24 + u32 v24 ), v̇2 = −u2 v2 (−a2 + 2a2 v2 − a2 v22 + a2 u22 v22 − u2 v23 − 2a2 u22 v23 + a2 u22 v24 ). Eliminating the common factor u2 between u̇2 and v̇2 rescaling the time we obtain the system (15) u̇2 = u2 (−a2 + 2a2 v2 − a2 v22 + a2 u22 v22 − u2 v23 − 2a2 u22 v23 − u32 v23 + a2 u22 v24 +u32 v24 ), v̇2 = −v2 (−a2 + 2a2 v2 − a2 v22 + a2 u22 v22 − u2 v23 − 2a2 u22 v23 + a2 u22 v24 ). The singular points of system (15) on the straight line u2 = 0 are (0, 0) and (0, 1). The eigenvalues of the linear part of system (15) evaluated at the singular point (0, 0) are −a2 and a2 , and so it is a hyperbolic saddle. The eigenvalues of the linear part of system (15) evaluated at the singular point (0, 1) are both 0 but the linear part of system (15) evaluated at (0, 1) is not linearly zero, so it is a nilpotent singular point. Before studying it we translate the point (0, 1) to the origin by doing the change (u2 , v2 ) = (u3 , 1 + v3 ). Therefore system (15) in the variables (u3 , v3 ) becomes u̇3 = −u23 − 3u23 v3 + u43 v3 − a2 u3 v32 − 3u23 v32 + a2 u33 v32 + 3u43 v32 − u23 v33 + 2a2 u33 v33 +3u43 v33 + a2 u33 v34 + u43 v34 , v̇3 = u3 + 4u3 v3 + a2 v32 + 6u3 v32 − a2 u23 v32 + a2 v33 + 4u3 v33 − 3a2 u23 v33 + u3 v34 −3a2 u23 v34 − a2 u23 v35 . Applying [4, Theorem 3.5] we get that it is a nilpotent saddle. Therefore the local phase portrait near the straight line u2 = 0 for system (15) is topologically equivalent to the one of Figure 3(a). Undoing the rescaling dt1 = u2 dt we get the local phase portrait of Figure 3(b). Going back through the changes of variables from the phase portrait of Figure 3(b), we obtain the local phase portrait at the origin of system (14) which is topologically equivalent to the one of Figure 3(c). Again going back through the changes of variables from the phase portrait of Figure 3(c), we obtain the local phase portrait around the origin of system (13) which is topologically equivalent to the one of Figure 3(d). Hence the origin of the local chart U1 has a nilpotent 12 J. LLIBRE AND C. VALLS saddle having two separatrices (one stable and one unstable) on the infinite circle and two separatrices (one stable and one unstable) in v < 0. (a) (b) (c) (d) Figure 3. Figures of the blow up of the singular point located at the origin of the local chart U1 of system (13). 4.2.2. The local chart U2 . System (1) in the chart U2 writes (16) u̇ = a2 u3 − a2 uv 2 + v 3 + u2 v 3 , v̇ = v(a2 u2 − a2 v 2 + uv 3 ). The origin is an infinite singular point in this chart. The linear part of system (16) at the origin is identically zero, so in order to determine its local phase portrait we must do blow ups. The characteristic directions at the origin of U2 are given by the real linear factors of v 4 = 0. Since u = 0 is not a characteristic direction, we do the vertical blow up (u, v) → (u1 , u1 v1 ), and system (16) writes in the new variables u̇1 = −u31 (−a2 + a2 v12 − v13 − u21 v13 ), v̇1 = −u21 v14 . Eliminating the common factor u21 between u̇1 and v̇1 rescaling the time we obtain the system (17) u̇1 = −u1 (−a2 + a2 v12 − v13 − u21 v13 ), v̇1 = −v14 . The unique singular point of system (17) on the straight line u1 = 0 is (0, 0). The eigenvalues of the linear part of system (17) at the singular point (0, 0) are 0 and a2 and so it is semi-hyperbolic. Applying [4, Theorem 2.19] we obtain that it is a saddle-node. Therefore the local phase portrait near the straight line u1 = 0 for system (17) is topologically equivalent to the one of Figure 4(a). Undoing the rescaling dt1 = u21 dt we get the phase portrait of Figure 4(b). Going back through the changes of variables from the phase portrait of Figure 4(b), we obtain the local phase portrait at the origin of system (16) which is topologically equivalent to the one of Figure 4(c). Hence the origin of the local chart U2 is topologically equivalent to an unstable node. 4.3. The infinite singular points of the local charts U1 and U2 when a2 a1 ̸= 0. We shall study the infinite singular points of the differential system (1) in the local charts U1 and U2 . We note that in this case doing the change of variables (x, y) → (−x, −y) we can assume, without loss of generality, that a2 a1 < 0. GLOBAL PHASE PORTRAITS OF THE GENERALIZED VAN DER POL SYSTEMS 13 (a) (b) (c) Figure 4. Figures of the blow up of the singular point located at the origin of the local chart U2 of system (16). 4.3.1. The local chart U1 . System (1) in the chart U1 writes (18) u̇ = −a2 u2 − a1 uv + a2 u2 v 2 − v 3 + a1 uv 3 − u2 v 3 , v̇ = −uv 4 . The unique infinite singular point in this chart is the origin. The linear part of system (18) at the origin is identically zero, so in order to determine its local phase portrait we must do blow ups. The characteristic directions at the origin are the real linear factors of uv(a2 u + a1 v) = 0. Since u = 0 is a characteristic direction before doing a vertical blow up we translate the direction u = 0 to u = v doing the change of variables (u, v) = (u1 − v1 , v1 ). Then system (18) becomes u̇1 = −a2 u21 − a1 u1 v1 + 2a2 u1 v1 + a1 v12 − a2 v12 + a2 u21 v12 − v13 + a1 u1 v13 (19) −2a2 u1 v13 − u21 v13 − a1 v14 + a2 v14 + u1 v14 , v̇1 = −u1 v14 + v15 . Now we do the vertical blow up (u1 , v1 ) → (u2 , u2 v2 ) and system (19) writes in the new variables u̇2 = u22 (−a2 − a1 v2 + 2a2 v2 + a1 v22 − a2 v22 + a2 u22 v22 − u2 v23 + a1 u22 v23 − 2a2 u22 v23 −u32 v23 − a1 u22 v24 + a2 u22 v24 + u32 v24 ), v̇2 = −u2 v2 (−a2 − a1 v2 + 2a2 v2 + a1 v22 − a2 v22 + a2 u22 v22 − u2 v23 + a1 u22 v23 − 2a2 u22 v23 −a1 u22 v24 + a2 u22 v24 ). Eliminating the common factor u2 between u̇2 and v̇2 rescaling the time we obtain the system u̇2 = u2 (−a2 − a1 v2 + 2a2 v2 + a1 v22 − a2 v22 + a2 u22 v22 − u2 v23 + a1 u22 v23 −2a2 u22 v23 − u32 v23 − a1 u22 v24 + a2 u22 v24 + u32 v24 ), (20) v̇2 = −v2 (−a2 − a1 v2 + 2a2 v2 + a1 v22 − a2 v22 + a2 u22 v22 − u2 v23 + a1 u22 v23 −2a2 u22 v23 − a1 u22 v24 + a2 u22 v24 ). The singular points of system (20) on the straight line u2 = 0 are (0, 0), (0, 1) and (0, a2 /(a2 − a1 )). The eigenvalues of the linear part of system (20) evaluated at the singular point (0, 0) are −a2 and a2 , and so it is a hyperbolic saddle. The eigenvalues of the linear part of system (20) evaluated at the singular point (0, 1) are 0 and −a1 and so it is semi-hyperbolic. Before studying it we translate the 14 J. LLIBRE AND C. VALLS point (0, 1) to the origin by doing the change (u2 , v2 ) = (u3 , 1 + v3 ). Therefore system (20) in the variables (u3 , v3 ) becomes u̇3 = −u23 + a1 u3 v3 − 3u23 v3 − a1 u33 v3 + u43 v3 + a1 u3 v32 − a2 u3 v32 − 3u23 v32 − 3a1 u33 v32 +a2 u33 v32 + 3u43 v32 − u23 v33 − 3a1 u33 v33 + 2a2 u33 v33 + 3u43 v33 − a1 u33 v34 + a2 u33 v34 + u43 v34 , v̇3 = u3 − a1 v3 + 4u3 v3 + a1 u23 v3 − 2a1 v32 + a2 v32 + 6u3 v32 + 4a1 u23 v32 − a2 u23 v32 − a1 v33 +a2 v33 + 4u3 v33 + 6a1 u23 v33 − 3a2 u23 v33 + u3 v34 + 4a1 u23 v34 − 3a2 u23 v34 + a1 u23 v35 −a2 u23 v35 . Now applying [4, Theorem 2.19] we get that it is a semi-hyperbolic saddle. The eigenvalues of the linear part of system (20) evaluated at the singular point (0, a2 /(a2 − a1 )) are 0 and a1 a2 /(a2 − a1 ), and so it is semi-hyperbolic. Before studying it we translate that point to the origin by doing the change (u2 , v2 ) = (u3 , a2a−a 2 1 + v3 ). Therefore system (20) in the variables (u3 , v3 ) after multiplied by (a2 − a1 )4 doing a convenient rescaling, we obtain u̇3 = u3 a1 a32 u3 − a42 u3 + a1 a32 u33 − a51 v3 + 4a41 a2 v3 − 6a31 a22 v3 + 4a21 a32 v3 − a1 a42 v3 −3a21 a22 u3 v3 + 6a1 a32 u3 v3 − 3a42 u3 v3 + a31 a22 u23 v3 − 2a21 a32 u23 v3 + a1 a42 u23 v3 −3a21 a22 u33 v3 + 2a1 a32 u33 v3 + a42 u33 v3 + a51 v32 − 5a41 a2 v32 + 10a31 a22 v32 − 10a21 a32 v32 +5a1 a42 v32 − a52 v32 + 3a31 a2 u3 v32 − 9a21 a22 u3 v32 + 9a1 a32 u3 v32 − 3a42 u3 v32 − 2a41 a2 u23 v32 +5a31 a22 u23 v32 − 3a21 a32 u23 v32 − a1 a42 u23 v32 + a52 u23 v32 + 3a31 a2 u33 v32 − 3a21 a22 u33 v32 −3a1 a32 u33 v32 + 3a42 u33 v32 − a41 u3 v33 + 4a31 a2 u3 v33 − 6a21 a22 u3 v33 + 4a1 a32 u3 v33 − a42 u3 v33 +a51 u23 v33 − 2a41 a2 u23 v33 − 2a31 a22 u23 v33 + 8a21 a32 u23 v33 − 7a1 a42 u23 v33 + 2a52 u23 v33 − a41 u33 v33 +6a21 a22 u33 v33 − 8a1 a32 u33 v33 + 3a42 u33 v33 − a51 u23 v34 + 5a41 a2 u23 v34 − 10a31 a22 u23 v34 +10a21 a32 u23 v34 − 5a1 a42 u23 v34 + a52 u23 v34 + a41 u33 v34 − 4a31 a2 u33 v34 + 6a21 a22 u33 v34  −4a1 a32 u33 v34 + a42 u33 v34 , v̇3 = (a2 + (a2 − a1 )v3 ) − a32 u3 + a41 v3 − 3a31 a2 v3 + 3a21 a22 v3 − a1 a32 v3 + 3a1 a22 u3 v3 −3a32 u3 v3 − a21 a22 u23 v3 + a1 a32 u23 v3 − a41 v32 + 4a31 a2 v32 − 6a21 a22 v32 + 4a1 a32 v32 − a42 v32 −3a21 a2 u3 v32 + 6a1 a22 u3 v32 − 3a32 u3 v32 + 2a31 a2 u23 v32 − 3a21 a22 u23 v32 + a42 u23 v32 + a31 u3 v33 −3a21 a2 u3 v33 + 3a1 a22 u3 v33 − a32 u3 v33 − a41 u23 v33 + a31 a2 u23 v33 + 3a21 a22 u23 v33 − 5a1 a32 u23 v33  +2a42 u23 v33 + a41 u23 v34 − 4a31 a2 u23 v34 + 6a21 a22 u23 v34 − 4a1 a32 u23 v34 + a42 u23 v34 . Now applying [4, Theorem 2.19] we get that it is a saddle-node. In this case the local phase portrait near the straight line u2 = 0 for system (20) is topologically equivalent to the one of Figure 5(a). Undoing the rescaling dt1 = u2 dt we get the phase portrait of Figure 5(b). Going back through the changes of variables from the phase portrait of Figure 5(b), we obtain the local phase portrait of system (19) in Figure 5(c). Finally, going back through the changes of variables from the local phase of Figure 5(c), we obtain the local phase portrait at the origin of system (18) which is topologically equivalent to the one of Figure 5(d). Hence the origin of the GLOBAL PHASE PORTRAITS OF THE GENERALIZED VAN DER POL SYSTEMS 15 local chart U1 is formed by four hyperbolic sectors (two stable and two unstable) separated by one parabolic sector. (a) (b) (c) (d) Figure 5. Figures of the blow up of the singular point located at the origin of the local chart U1 of system (18). 4.3.2. The local chart U2 . System (1) in the chart U2 writes (21) u̇ = a2 u3 +a1 u3 v−a2 uv 2 +v 3 −a1 uv 3 +u2 v 3 , v̇ = a2 u2 v+a1 u2 v 2 −a2 v 3 −a1 v 4 +uv 4 . The origin is an infinite singular point in this chart. The linear part of system (21) at the origin is identically zero, so in order to determine its local phase portrait we must do blow ups. The characteristic directions at the origin of U2 are given by the real linear factors of v 4 = 0. Since u = 0 is not a characteristic direction, we do the vertical blow up (u, v) → (u1 , u1 v1 ) and system (21) writes in the new variables u̇1 = −u31 (−a2 − a1 u1 v1 + a2 v12 − v13 + a1 u1 v13 − u21 v13 ), v̇1 = −u21 v14 . Eliminating the common factor u21 between u̇1 and v̇1 rescaling the time we obtain the system (22) u̇1 = −u1 (−a2 − a1 u1 v1 + a2 v12 − v13 + a1 u1 v13 − u21 v13 ), v̇1 = −v14 The unique singular point of system (22) on the straight line u1 = 0 is (0, 0). The eigenvalues of the linear part of system (22) at the singular point (0, 0) are 0 and a2 and so it is semi-hyperbolic. Applying [4, Theorem 2.19] we obtain that it is a saddle-node. Going back through the changes of variables and reascaling of the time we get that the origin of the local chart U2 is topologically equivalent to the one provided in Figure 4(c) when a2 > 0 and Figure 4(c) reversing the orientation of the orbits when a2 < 0. Proof of Theorem 2. When a2 = 0 and a1 > 0 the two stable separatrices of the origins of U1 and V1 are at infinity, and their unstable separatrices go to the finite part of the Poincaré disc. Moreover the origins of the local charts U2 and V2 are unstable nodes. Since the origin is an unstable focus if a1 ∈ (0, 2) and an unstable node if a1 ≥ 2 the Poincaré Bendixson Theorem forces the existence of 2 √ f (x) = a1 (x − 1), g(x) = x, at least one limit cycle. In view of Theorem 6 with 3 F (x) = a1 (x /3 − x), a = −1, b = 1 and c = 3, we get that the limit cycle is unique. Combining all this information we obtain that the global phase portrait of 16 J. LLIBRE AND C. VALLS system (1) when a2 = 0 is topologically equivalent to the one of Figure 1(a) and so statement (i) is proved. When a1 = 0 and a2 > 0 taking into account the local phase portraits at the origin of coordinates (which is a center), at the origins of U1 and V1 that are formed by the four hyperbolic sectors of the nilpotent saddle of Figure 3(d), the unstable node at the origin of U2 , the stable node at the origin of V2 , and that the global phase portrait is symmetric with respect to the x-axis (because the differential system (1) is invariant under the change (x, y, t) → (x, −y, −t)) we get that the phase portrait in the Poincaré disc is topologically equivalent to the one of Figure 1(b). Hence statement (ii) is proved. Finally, when a1 a2 ̸= 0 and a1 a2 < 0 with a2 > 0, taking into account the local phase portraits at the origin of coordinates (which is a either a stable focus or a stable node), at the origins of U1 and V1 that are formed by the four hyperbolic sectors separated by parabolic sectors of Figure 5(d), the unstable node at the origin of U2 and the stable node at the origin of V2 we get that a possible phase portrait in the Poincaré disc is topologically equivalent to the one of Figure 1(c) and we see that there are numerical values such that this phase portrait is realizable. Note that in this phase portrait there is a limit cycle. Hence the theorem is proved. In the case of a2 < 0 we get the same global phase portrait reverting the orientation of the orbits. There are topologically two more possible global phase portraits but numerical evidence seem to show that they are not realizable. □ Acknowledgements The first author is supported by the Agencia Estatal de Investigación grant PID2019-104658GB-I00 and the H2020 European Research Council grant MSCA- RISE-2017-777911. The second author is partially supported by FCT/Portugal through CAMGSD, IST-ID, projects UIDB/04459/2020 and UIDP/04459/2020. References [1] M.J. Álvarez, A. Ferragud and X. Jarque, A survey on the blow up technique, Int. J. Bifur. Chaos 21 (2011), 3103–3118. [2] A.A. Andronov, E.A. Leontovich, I.I. Gordon and A.G. Muier, Qualitative theory of second-order dynamic systems, John Wiley & Sons, New York, Toronto, 1973. [3] J. Burke, M. Desroches, A. Granados, T.J. Kaper, M. Krupa and T. Vo, From canards of folded singularities to torus canards in a forced van der Pol equation, J. Nonlinear Sci. 26 (2016), no. 2, 405–451. [4] F. Dumortier, J. Llibre and J.C. Artés, Qualitative theory of planar differential systems, UniversiText, Springer-Verlag, New York, 2006. [5] M. Hirsch, C. Pugh and M. Shub, Invariant Manifolds, Lecture Notes in Mathematics, Vol. 583. Springer-Verlag, 1977. [6] N. Levinson and O.K. Smith, A general equation for relaxation oscillations, Duke Math. J. 9 (1942), 382–403. [7] X. Ma, Q. Bi and L. Wang, Complex periodic bursting structures in the Rayleigh–van der Pol–Duffing oscillator J. Nonlinear Sci. 32 (2022), no. 2, Paper No. 25, 23 pp. [8] Z. Wang and H. Chen, A nonsmooth van der Pol-Duffing oscillator (I): The sum of indices of equilibria is −1, Discrete Contin. Dyn. Syst. Ser. B 27 (2022), no. 3, 1421–1446. [9] Z. Wang and H. Chen, A nonsmooth van der Pol-Duffing oscillator (II): The sum of indices of equilibria is 1, Discrete Contin. Dyn. Syst. Ser. B 27 (2022), no. 3, 1549–1589. GLOBAL PHASE PORTRAITS OF THE GENERALIZED VAN DER POL SYSTEMS 17 [10] Zhang Zhi–Fen et al, Qualitative Theory Differential Equations, Translation of Mathemat- ical Monographs, 101, American Mathematical Society, 1991. 1 Departament de Matemàtiques, Universitat Autònoma de Barcelona, 08193 Bel- laterra, Barcelona, Catalonia, Spain Email address: jllibre@mat.uab.cat 2Departamento de Matemática, Instituto Superior Técnico, Universidade de Lisboa, Av. Rovisco Pais 1049–001, Lisboa, Portugal Email address: cvalls@math.ist.utl.pt