This is a preprint of: “A family of periodic orbits for the extended Hamiltonian system of the Van der Pol oscillator”, Jean-Marc Ginoux, Jaume Llibre, J. Geom. Phys., vol. 183, 104705, 2023. DOI: [10.1016/j.geomphys.2022.104705] A FAMILY OF PERIODIC ORBITS FOR THE EXTENDED HAMILTONIAN SYSTEM OF THE VAN DER POL OSCILLATOR JEAN-MARC GINOUX1 AND JAUME LLIBRE2 Abstract. In this paper we show a new way of using the averag- ing theory for studying families of periodic orbits of a Hamiltonian system. We do this study computing a new family of periodic or- bits of the extension of the Van der Pol oscillator to a Hamiltonian system of two degrees of freedom. 1. Introduction and statement of the main results The classical Van der Pol oscillator is modeled by the second-order differential equation (1) ẍ − µ(1 − x2 )ẋ + x = 0, where x = x(t) is the position coordinate at time t, µ is a parameter, and the dot denotes derivative with respect to the time t, see [9, 10, 11]. Initially the differential equation (1) allowed to Van der Pol to ex- plain the stable oscillations observed in electrical circuits employing vacuum tubes. Later on this differential equation has been used for explaining different phenomena in biology, physics, . . . Thus in biology was utilised as a model for studying the action potentials of neurons, see for instance [3, 4, 7]. While in physics equation (1) was also used in phonation to model the right and left vocal fold oscillators (see [6]), and in seismology to model two plates in a geological fault (see [1]), . . . More recently in 2015, see equations (9) of [8], the Van der Pol oscillator was written in the Hamiltonian formalism by extending it to a four-dimensional autonomous differential as follows ẍ − µ(1 − x2 )ẋ + x = 0, (2) ÿ + µ(1 − x2 )ẏ + y = 0. 2010 Mathematics Subject Classification. 37J45. Key words and phrases. periodic orbit, Hamiltonian system, Van der Pol oscil- lator, averaging theory. 1 2 J.M. GINOUX, J. LLIBRE Observe that the motion of the Van der Pol oscillator (1) is not affected because the evolution of the variable x(t) is independent of the variable y(t). Defining py = ẋ and px = ẏ + µ(1 − x2 )y from (2) we obtain that this differential system of two differential equations of second order admits the following Hamiltonian formalism ẋ = py , ẏ = px − µ(1 − x2 )y, (3) p˙x = −y − µ2xypy , p˙y = −x + µ(1 − x2 )py , with the Hamiltonian (4) H = px py + xy − µ(1 − x2 )ypy . Of course px and py are the conjugate momenta of the variables x and y, respectively. In [2] it is shown that the Hamiltonian system (3) establishes a connection between the phase of the limit cycle of the Van der Pol system with the Hannay angle of the Hamiltonian system. In this paper our objective is to study analytically the periodic orbits of the Hamiltonian system (3) using the averaging theory. In general the periodic orbits of the Hamiltonian systems are studied numerically because its analytical study is in general difficult and sometimes im- possible. Here the technique used with the averaging theory can be extended to many other Hamiltonian systems. Our main result is the following. Theorem 1. For µ sufficiently small the following statements hold for the Hamiltonian system (3). (a) For each h ∈ R there exists a periodic orbit  γh (t) = x(t; µ, h), y(t; µ, h), px (t; µ, h), py (t; µ, h) such that when the parameter µ → 0 satisfies that γh (t) tend to the periodic orbit p ! p ! 4 − µ2 |h| ∗ 4 − µ2 2 cos t , cos φ + t , 2 2 2 p ! p !! |h| 4 − µ 2 4 − µ 2 − sin φ∗ + t , −2 sin t 2 2 2 of the Hamiltonian system (3) with µ = 0. 3 (b) The periodic orbit γh (t) live on the energy level H = h, and when h varies in R we get a family of periodic orbits. Theorem 1 is proved in the section 3. In section 2 we present the algorithm that we follow for studying the families of periodic orbits of the Hamiltonian systems. Since the differential system formed by the equations ẋ and p˙y is the first order differential system associated to the Van der Pol oscillator (1), and this differential system does not depend on the variables y and py , it follows immediately from Theorem 1 the following corollary. Corollary 2. For µ sufficiently small and for each h ∈ R the Van der Pol oscillator has a periodic orbit  γh (t) = x(t; µ, h), py (t; µ, h) such that when the parameter µ → 0 satisfies that γh (t) tend to the periodic orbit p ! p !! 4 − µ2 4 − µ2 2 cos t , −2 sin t , 2 2 solution of the Van der Pol oscillator with µ = 0. Of course, for a fixed value of µ sufficiently small the periodic orbit γh (t) of Corollary 2 is the classical limit cycle of the Van der Pol oscil- lator, because it is known that such a limit cycle is unique for a given value of µ. 2. The algorithm The algorithm that we will follow for studying the families of periodic orbits of a Hamiltonian system has the following steps. 1. We select an equilibrium point of the Hamiltonian system having at least one pair of complex conjugate eigenvalues and translate it to the origin of coordinates. 2. We do a change of variables which writes the matrix of the linear part of the Hamiltonian system at the origin of coordinates in its real Jordan normal form. Then, in general the obtained differential system have lost its Hamiltonian structure. 3. The coordinates of the differential system whose linear part at the origin of coordinates produce complex eigenvalues are changed to polar 4 J.M. GINOUX, J. LLIBRE coordinates. More precisely, if the real Jordan normal form of the linear part at the origin of coordinates has the block   α −β , β α associated to the coordinates x and y with α non identically zero. Then we do the change of variables (x, y) → (r, θ), where x = r cos θ and y = r sin θ. We shall take α as the small parameter µ necessary for applying the averaging theory, see the appendix. But if there are more than one block associated to complex eigenval- ues, for all the other blocks different from the first we do the following change of variables. Take one of these blocks, for instance   a −b , b a associated to the coordinates u and v and do the change of variables (u, v) → (R, φ), where x = R cos(φ + kθ), y = R sin(φ + kθ), and the constant k is chosen in such a way that the differential system in the new coordinates must satisfy φ̇ = O(α), with α sufficiently small. This step will become more clear when we apply the algorithm to the Hamiltonian system (3) in the next section. 4. Now we take the angle variable θ as the new independent variable. 5. Since generically the periodic orbits in the Hamiltonian systems appear in families of periodic orbits parametrized by h, being h the different values of the energy levels H = h of the Hamiltonian H, and the averaging theory only detects isolated orbits, we must apply the averaging theory in the energy levels H = h, eliminating one variable using the relation H = h. If the Hamiltonian H is the unique inde- pendent first integral of the Hamiltonian system the differential system obtained on the energy level H = h is into the normal form (14) of the appendix, ready for applying to it the averaging theory for studying its periodic orbits. 6. In case that the Hamiltonian system has additional independent first integrals Fk , also independent with the Hamiltonian H of the Hamiltonian system, we need to fix them, i.e. Fk = fk and work on these invariant subspaces in order that we can apply the averaging theory and to study the periodic orbits of the Hamiltonian system. 5 3. Proof of Theorem 1 Now we apply the algorithm described in section 2 to the Hamilton- ian system (3) for proving Theorem 1. We will indicate all the steps of the algortihm using its numbers. 1. It is easy to verify that the Hamiltonian system (3) has a unique equilibrium point localized p at the origin of coordinates, and that its four eigenvalues are (±µ ± µ2 − 4)/2. 2. Now we write the matrix   0 0 0 1  0 −µ 1 0     0 −1 0 0  −1 0 0 µ of linear part of the Hamiltonian system (3) into its real Jordan normal form  µ 1p  − 4 − µ2 0 0  p 2 2   2   4−µ µ   0 0   2 2   µ 1 p ,  2   0 0 − − 4 − µ   p 2 2   2  4−µ µ 0 0 − 2 2 doing the change of variables (x, y, px , py ) → (X, Y, Z, V ), where   1 p 2  1 p 2 (x, y, px , py ) = X, Z, µZ − 4 − µ V , µX − 4 − µ Y . 2 2 Then the Hamiltonian system (3) becomes p µ 4 − µ2 Ẋ = X − Y, 2 2 p 4 − µ2 µ µ2 Ẏ = X+ Y +p X 3 − µX 2 Y, 2 2 4−µ 2 (5) p 4 − µ2 µ Ż = − V − Z + µX 2 Z, 2 2 p µ 4 − µ2 3µ2 V̇ = − V + Z+p X 2 Z − 2µXY Z, 2 2 4 − µ2 6 J.M. GINOUX, J. LLIBRE and the Hamiltonian (4) becomes the first integral 1 p  p  H= µZ − 4 − µ2 V µX − 4 − µ2 Y − (6) 4 1  p  µ (1 − X 2 ) Z µX − 4 − µ2 Y + XZ. 2 of the differential system (5). 3. We do the change of variables (X, Y, Z, V ) → (r, θ, R, φ) where  (X, Y, Z, V ) = r cos θ, r sin θ, R cos(φ + kθ), R sin(φ + kθ) , and the differential system (5) writes in the new variables as ! ! 1 µ cos θ ṙ = µ r + sin θ cos2 θ p − sin θ r3 , 2 4−µ 2 p ! 4 − µ2 µ cos θ θ̇ = + µ cos3 θ p − sin θ r2 , 2 4 − µ2 1 1 6µ cos θ sin(φ + θ) Ṙ = µ − R + cos θ cos(φ + θ) p + (7) 2 2 4 − µ2 3 cos(φ + 2θ) − cos φ) r2 R) , ! 3µ cos(φ + θ) φ̇ = µ cos θ cos θ cos(φ + θ) p − sin(φ + θ) 4 − µ2 ! µ cos3 θ −p − 2 sin θ cos2 (φ + θ) + sin θ cos2 θ r2 , 4 − µ2 where we have take k = 1 in order that φ̇ = O(µ). The first integral H given in (6) in the variables (r, θ, R, φ) writes 1 H = rR cos θ cos(φ + θ) − µR cos(φ + θ) (1 − r2 cos2 θ) 2  p  µr cos θ − p4 − µ2 r sin θ (8) µr cos θ − 4 − µ2 r sin θ + 4  p  2 µR cos(φ + θ) − 4 − µ R sin(φ + θ) . 7 4. We take now the variable θ as the new independent variable and the differential system (7) becomes (9) µ  r′ = r − 2r3 sin2 θ cos2 θ + O(µ2 ), 2 µ R′ = − (1 + r2 cos θ cos(φ + θ)(cos φ − 3 cos(φ + 2θ))) R + O(µ2 ), 2 2 µr φ′ = cos θ(sin(2φ + θ) − 3 sin(2φ + 3θ) − 3 sin θ + sin(3θ)) + O(µ2 ), 4 where the prime denotes derivative with respect to the variable θ. 5. From (8) and H = h we obtain that   h (10) φ = arccos + O(µ). rR Substituting this φ given in (10) we obtain the following expression for the differential system (9) restricted to the energy level H = h (11) r  r′ = µ − r3 sin2 θ cos2 θ + O(µ2 ) 2 = µF11 (θ, r, R) + O(µ2 ),       ′ µ h R = −hr cos θ cos arccos +θ −R + 2 rR         2 h h 3r R cos θ cos arccos + θ cos arccos + 2θ rR rR +O(µ2 ) = µF12 (θ, r, R) + O(µ2 ). The differential system (11) is already written in the normal form (14) of the appendix for applying the averaging theory. Here us- ing the notation of the appendix we have that F = (F11 , F12 ), x = (r, R) and T = 2π then we compute the averaged function f (x) = (f11 (r, R), f12 (r, R)) defined in (16) of the appendix, and we get   1 r3 f11 (r, R) = r− , 2 4 (12) 1 2 h2 f12 (r, R) = (3r − 4) R − . 8 4R The real solutions of the system f11 (r, R) = 0 and f12 (r, R) = 0 with r and R positive are r = 2 and R = |h|/2. Moreover, on this solution 8 J.M. GINOUX, J. LLIBRE the Jacobian (17) for our differential system (12) is (13)   1 (4 − 3r2 ) 0  8  det    2    = −8π 2 ̸= 0. 3rR 1 2h + 3r2 − 4 4 8 R2 (r,R)=(2,|h|/2) From Theorem 3 of the appendix we get for µ sufficiently small that for each h ∈ R we have the periodic solution  γh1 (θ; µ, h) = r(θ; µ, h), R(θ; µ, h) of the differential system (12) such that γh1 (θ; µ, h) → (2, |h|/2), when µ → 0. The periodic orbit γh1 (θ; µ, h) provides the following periodic orbit of the differential system (9)  γh2 (θ; µ, h) = r(θ; µ, h), R(θ; µ, h), φ(θ; µ, h) , such that γh2 (θ; µ, h) → (2, |h|/2, φ∗ ), when µ → 0, where φ∗ = 0 if h > 0 and φ∗ = π if h < 0. The periodic orbit γh2 (θ; µ, h) provides the following periodic orbit of the differential system (7)  γh3 (t; µ, h) = r(t; µ, h), θ(t; µ, h), R(t; µ, h), φ(t; µ, h) ,  p  such that γh3 (t; µ, h) → 2, 4 − µ2 t/2, |h|/2, φ∗ , when µ → 0. The periodic orbit γh3 (t; µ, h) provides the following periodic orbit of the differential system (5)  γh4 (t; µ, h) = X(t; µ, h), Y (t; µ, h), Z(t; µ, h), V (t; µ, h) , such that p ! p ! 4 − µ2 4 − µ2 γh4 (t; µ, h) → 2 cos t , 2 sin t , 2 2 p ! p !! |h| ∗ 4 − µ2 |h| ∗ 4 − µ2 cos φ + t , sin φ + t , 2 2 2 2 when µ → 0. The periodic orbit γh4 (t; µ, h) provides the following periodic orbit of the differential system (3)  γh (t; µ, h) = x(t; µ, h), y(t; µ, h), px (t; µ, h), py (t; µ, h) , 9 such that p ! p ! 4 − µ2 |h| 4 − µ 2 γh (t; µ, h) → 2 cos t , cos φ∗ + t , 2 2 2 p ! p !! |h| 4 − µ 2 4 − µ 2 − sin φ∗ + t , −2 sin t , 2 2 2 when µ → 0. This completes the proof of Theorem 1. Appendix Here we summarize the basic results from averaging theory for com- puting periodic orbits that we shall use for proving Theorem 1, for a proof see Theorem 11.5 of [12], or Theorem 1.2.1 of [5]. First we consider the initial value problem (14) ẋ = µF (θ, x) + µ2 R(θ, x, µ), x(0) = x0 , with x ∈ U where U is an open subset of Rn , and θ ≥ 0. We assume that F (θ, x) and R(θ, x, µ) are periodic functions in θ of period T . Of course the dot in (14) denotes derivative with respect to the variable θ. Second we consider in the open set U the following initial value problem for the averaged differential system (15) ẏ = µf (y), y(0) = x0 , where f (x) is the averaged function of the function F (θ, x), i.e. Z 1 T (16) f (x) = F (θ, x)dθ. T 0 Theorem 3. For the two initial value problems (14) and (15) we as- sume: (i) the function F is C 2 and bounded by a constant which does not depend on µ in [0, ∞) × U for all µ ∈ (0, µ0 ]. (ii) the functions F and R are periodic in θ of period T , and T is independent of µ. (iii) the solution y(θ) of (15) belongs to U for all θ ∈ [0, 1/µ]. Then for each equilibrium point p of the averaged system (15) satisfying   ∂f (17) det ̸= 0, ∂y y=p 10 J.M. GINOUX, J. LLIBRE there is a periodic solution x(θ, µ) of period T for the system (14) verifying x(θ, µ) → p as µ → 0. Declaration of Interest The results of this paper are new and original. They do not have any conflict of interest. Acknowledgements This work is partially supported by the Agencia Estatal de Inves- tigación grant PID2019-104658GB-I00, and the H2020 European Re- search Council grant MSCA-RISE-2017-777911. References [1] J.H.E. Cartwright, V.M. Eguı́luz, E. Hernández-Garcı́a and O. Piro, Dynamics of elastic excitable media, Int. J. Bifurcation and Chaos 9(11) (1999), 2197- 2202. [2] R. Chattopadhyay, T. Shah, S. Chakraborty, Finding the Hannay angle in dissipative oscillatory systems via conservative perturbation theory Physical Review E 97(6) (2018), 062209. [3] R. FitzHugh, Impulses and physiological states in theoretical models of nerve membrane, Biophysical Journal 1 (6) (1961), 445-466. [4] J.M. Ginoux, History of Nonlinear Oscillations Theory, Archimede, New Stud- ies in the History and Philosophy of Science and Technology, vol. 49, Springer, New York, 2017. [5] J. Llibre, R. Moeckel and C. Simó, Central configurations, periodic orbits and Hamiltonian systems, Advances Courses in Mathrematics, CRM Barcelona, Birhauser, 2015. [6] J. Lucero and J. Schoentgen, Modeling vocal fold asymmetries with coupled van der Pol oscillators, Proceedings of Meetings on Acoustics, Vol. 19, 2013, 060165. [7] J. Nagumo, S. Arimoto and S. Yoshizawa, An active pulse transmission line simulating nerve axon, Proceedings of the IRE, Institute of Electrical and Electronics Engineers (IEEE) 50(10) (1962), 2061-2070. [8] T. Shah, R. Chattopadhyay, K. Vaidya and S. Chakraborty, Conservative perturbation theory for non-conservative systems, Phys. Rev. E 92, (2015), 062927. [9] B. van der Pol, A theory of the amplitude of free and forced triode vibrations, Radio Review 1 (1820), 701-710. [10] B. van der Pol, On relaxation-oscillations, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 2(11) (1926), 978-992. [11] B. van der Pol and J. van der Mark, Frequency demultiplication, Nature 120(3019) (1927), 363-364. 11 [12] F. Verhulst, Nonlinear Differential Equations and Dynamical Systems, Univer- sitext, Springer, 1991. 1 Aix Marseille Univ, Université de Toulon, CNRS, CPT, Marseille, France Email address: ginoux@univ-tln.fr 2 Departament de Matemàtiques, Universitat Autònoma de Barcelona, 08193 Bellaterra, Barcelona, Spain Email address: jllibre@mat.uab.cat