Global Dynamics and Bifurcation of Periodic Orbits in a Modified Nosé-Hoover Oscillator

We perform a global dynamical analysis of a modified Nosé-Hoover oscillator, obtained as the perturbation of an integrable differential system. Using this new approach for studying such an oscillator, in the integrable cases, we give a complete description of the solutions in the phase space, including the dynamics at infinity via the Poincaré compactification. Then using the averaging theory, we prove analytically the existence of a linearly stable periodic orbit which bifurcates from one of the infinite periodic orbits which exist in the integrable cases. Moreover, by a detailed numerical study, we show the existence of nested invariant tori around the bifurcating periodic orbit. Finally, starting with the integrable cases and increasing the parameter values, we show that chaotic dynamics may occur, due to the break of such an invariant tori, leading to the creation of chaotic seas surrounding regular regions in the phase space.


Introduction
Since its appearance in 1984, the now well-known Nosé-Hoover oscillator, given by the system of ordinary differential equations [11,12]: was widely studied. Obtained from the propositions made by Nosé in the seminal paper [11] as new paradigms in the study of molecular dynamics, the apparently simple system (1) was formulated and further studied by Hoover and collaborators in [4,13], showing its very rich dynamical behavior, as the existence of several types of periodic orbits, nested invariant tori, and even chaotic behavior. In synthesis, system (1) models the one-dimensional harmonic oscillator obtained using Nosé's canonical equations of motion. In this way, the time-dependent variables x, y, and z represent, respectively, the momentum, the position coordinate, and the friction coefficient of a particle. As the role of the friction coefficient z is to maintain the average temperature equal to 1, there is a control parameter, denoted by α. System (1) and some of its generalizations were studied by several authors, from mathematical and physical points of view (see for instance [3, 4, 11, 13-15, 17, 18, 20] and references therein). For details on the physics background for obtaining system (1), see [11] and the nice introduction and motivation section in [13]. More recently, in [15], the authors presented a brief but interesting review about the study of Nosé-Hoover oscillator throughout the years.
Observe that for α > 0 system (1) has no equilibrium points. Then, the standard methods often used in dynamical system theory, as the determination of equilibrium points, the study of their stability and bifurcations, seeking for periodic orbits bifurcating from them, or the determination of their stable and unstable manifolds and their possible transversal intersections, which may lead to chaotic dynamics, cannot be applied for studying the dynamics of system (1). Hence, other techniques must be employed in order to study it.
In [18], the authors studied existence of periodic orbits of system (1). In that paper, they said that there are two obvious approaches to the theory of Nosé-Hoover equations: one can start from the study of its trajectories for α small or, alternatively, one can study the trajectories which pass close to infinity. Using the second possible approach, in [18], the authors showed the existence of several types of periodic orbits of system (1) for 0 < α ≤ 1, which bifurcate from orbits heteroclinic to equilibrium points at infinity. We shall see ahead that with this approach at least one relevant periodic orbit, which does not bifurcate from infinity, was missed.
In this paper, we propose a third alternative to study system (1): we consider the Nosé-Hoover oscillator as a perturbation of an integrable differential system, which enables us to perform a global dynamical analysis of it. More precisely, we study the following differential system:ẋ where x, y, z ∈ R 3 are the state variables, a and b are real parameters, and the dot denotes derivative with respect to the independent variable t. System (2) was obtained from the Nosé-Hoover oscillator (1) by means of the rescaling z → aZ of the z variable. Doing this rescaling and then considering b = α/a and calling Z = z again, we obtain system (2), which for a = 0 is topologically equivalent to system (1).
The apparent artificial change of variable considered above to obtain system (2) from system (1) has a practical advantage: the obtained two-parameter system is completely integrable for a = 0 or b = 0. Although system (2) is not topologically equivalent to system (1) for a = 0, we can start the analysis of system (2) for a = 0 (or b = 0) and then study what happens for a, b = 0 small, that is, we consider perturbations of the integrable cases, which are equivalent to system (1). With this new approach, we intend to make a contribution to understanding the complicated dynamical behavior of Nosé-Hoover oscillator (1). In particular, using the averaging method, we analytically prove the existence of a linearly stable periodic orbit bifurcating from one of the infinite periodic orbits which exist in the integrable cases. This periodic orbit, which was missed in [18], is a fundamental dynamical element for the existence of nested invariant tori of system (2) (and, consequently, of Nosé-Hoover oscillator (1)), because the nested invariant tori exist exactly around this periodic orbit, as we shall see ahead. We also investigate numerically the occurrence of chaotic behavior for system (2), by increasing the values of parameters from a = 0 and b = 0.
The rest of this paper is organized as follows. In Section 2, we study system (2) in the integrable cases, that is, when a = 0 and b = 0. In Section 3, we study the global dynamics of system (2), including the dynamics at infinity via the Poincaré compactification for a polynomial vector field in R 3 . In Section 4, using the averaging method, we prove the existence of a periodic orbit bifurcating from the integrable systems. In Section 5, we perform a numerical analysis of system (2) for a and b different from zero, showing the existence of invariant tori around the bifurcated periodic orbit and the occurrence of chaotic behavior, due to the destruction of these invariant tori, as the parameter values are varied. Finally, in Section 6, we present some concluding remarks and comments.

Dynamics of System (2) in the Integrable Cases
Differential system (2) is completely integrable for a = b = 0; a = 0 and b = 0; a = 0 and b = 0. For a = b = 0, the dynamics of system (2) is trivial; the phase space is foliated by the invariant planes z = c, with c ∈ R. In this case, the z-axis is filled by the equilibrium points (0, 0, z), which are centers on the invariant planes z = c, as it is shown in Fig. 1.
For a = 0 and b = 0, system (2) has two independent first integrals, given by: and with x = 0. For b = 0 and a = 0, system (2) has the first integral H 2 (x, y, z) = z and with x = 0 and z = ±2/a. Due to the complete integrability of system (2) in these two cases, it is possible to determine completely the global dynamics of its solutions, which is described in the following theorems. Proof If a = 0, b > 0 and x = y = 0, then system (2) reduces to: Hence, the z-axis is invariant under the flow of system (2) and the dynamics on it is governed by the equationż = −b < 0, from which follows statement (a) of Theorem 1. For a = 0, system (2) has the independent first integrals H 1 and F 1 given by Eq. 3 and Eq. 4, respectively. Hence, the orbits are contained in the intersection {H 1 = c} ∩ {F 1 = k}, with c, k ∈ R. For c = 2, we have from Eq. 4 that the solutions are contained in the surfaces 2z − bxy = k, k ∈ R. Now, using Eq. 3, we have x = ± 2 − y 2 , from which follows that which gives a closed curve for each value of k ∈ R (see Fig. 2, which was obtained considering b = 1 and k = 0). Thus, the invariant cylinder x 2 + y 2 = 2 is filled by periodic orbits of system (2), as shown in Fig. 3 right. This proves statement (b) (i) of Theorem 1 If x 2 + y 2 = c with 0 < c = 2, then the orbits are contained in the intersection of the surfaces x 2 + y 2 = c and F 1 (x, y, z) = k, determined by the curves which, for each k ∈ R, gives an upward spiraling curve for c > 2 and a downward spiraling curve for 0 < c < 2, on the invariant cylinder x 2 + y 2 = c, see Fig. 3 left. Therefore, statements (b) (ii) and (iii) of Theorem 1 are proved.

Theorem 2
Consider system (2) with b = 0 and a > 0. In this case, the following statements hold.
(a) The z-axis is filled by equilibrium points of system (2).  then (0, 0, c) is a stable node; and if c < −(2/a), then it is an unstable node.
If b = 0 and a < 0 the same results hold, reversing appropriately the stability of the equilibrium points (0, 0, c).
Proof If b = 0, then the z-axis is filled by equilibrium points of system (2). Moreoverż = 0 implies that the phase space R 3 is foliated by the invariant planes z = c, with c ∈ R. In each invariant plane, system (2) reduces to the following linear differential system: where c ∈ R and a > 0. The eigenvalues of the linear part of system (6) at the origin are: from which follows trivially the statements of Theorem 2. The same analysis can be made for a < 0. The dynamics of system (2) in the integrable cases are of fundamental importance to describe the complicated dynamical elements that appear in the phase space of this system for a and b different from zero, as periodic orbits, nested invariant tori and chaotic behavior, as we shall see ahead.

Dynamics of System (2) at Infinity
In order to study the global dynamics of system (2), including the unbounded solutions, this polynomial differential system can be extended by a change of coordinates into an analytic differential system defined on a closed ball of radius 1 (the Poincaré ball), whose interior is diffeomorphic to R 3 and its boundary, the 2-dimensional sphere S 2 = {(x, y, z) : , plays the role of the infinity of R 3 . In this way, it is possible to lead the unbounded phase space R 3 of system (2) to a compact phase space given by the closed ball of radius 1. One of the known techniques for making such an extension is the Poincaré compactification for polynomial vector fields, which is described for instance in [1], and some applications of it can be found in [6,7,9]. The Poincaré compactification enables us to study the dynamics of system (2) near and at infinity and to analyze how the unbounded solutions come from or go to infinity.
When we perform the Poincaré compactification of system (2), we obtain six polynomial vector fields defined on the local charts U i and V i , i = 1, 2, 3, with coordinates (z 1 , z 2 , z 3 ) which cover the sphere as a differential manifold. All the points on the invariant sphere (at infinity) in the coordinates of any chart U i and V i are obtained simply by making z 3 = 0 in the transformed systems. The points in the interior of the Poincaré ball, which is diffeomorphic to R 3 , are given in the local charts U i by z 3 > 0 and in the local charts V i by z 3 < 0.  Fig. 4 for an illustration of the sphere at infinity S 2 and for the orientation of the local charts U i and V i . We use this compactification technique to study the orbits of system (2) near and at infinity. The expression of system (2) in the local chart U 1 , after making z 3 = 0 iṡ The z 2 -axis is invariant under the flow of system (7). If ab < 0, then system (7) has the equilibrium points: and the straight lines z 2 = ± √ −b/a are invariant under the flow of system (7). Calculating the eigenvalues of the linear part of system (7) at the equilibrium points P ± , we obtain that, for b < 0 < a, P + is an unstable node and P − is a stable node, as shown in Fig. 5a. In the case a < 0 < b, these nodes change their stability. Note that if ab > 0, system (7) has no equilibrium points; hence, the phase portrait is as shown in Fig. 5b. Now, we consider the compactification in the chart U 2 . After doing z 3 = 0, we get the following differential system at infinity: System (8) has a straight line filled by equilibrium points in the z 2 -axis. If b < 0 < a, then the equilibrium points (0, z 2 ) are normally hyperbolic, being stable for z 2 > 0 and unstable for z 2 < 0. Rescaling the system by z 1 = 0, we have that (0, 0) is of saddle type (see Fig. 6a). For ab > 0, the solutions are contained in the invariant ellipses: with k > 0 constant. These ellipses are formed by heteroclinic orbits to the equilibrium points in the z 2 -axis, as shown in Fig. 6b, considering a > 0 and b > 0. The cases a < 0 < b and a < 0 and b < 0 have the same phase portraits, reversing the orientation of the orbits. In the local chart U 3 , after doing z 3 = 0, we obtain the following system at infinity: Note that the z 2 -axis is filled by equilibrium points of system (9) and the z 1 -axis is invariant under the flow of this system. Moreover, for ab < 0, system (9) has the equilibrium points: while, for ab > 0, there are no equilibrium points outside the z 2 -axis. Studying the eigenvalues of the Jacobian matrix at these equilibrium points and rescaling system (9) by z 1 = 0, we obtain the phase portraits shown in Fig. 7, (a) for b < 0 < a and (b) for a > 0 and b > 0. In the case a < 0 < b and a < 0 and b < 0, the phase portraits are the same reversing the orientation of the orbits. The flow in the local charts V i , for i = 1, 2, 3, restricted to z 3 = 0 is the same as the flow in the respective local charts U i , reversing the time, because the compactified vector field in V i is the same as the vector field in U i multiplied by −1. Hence, the phase portraits are the same, reversing the orientation of the orbits.
Considering the analysis of system (2) at infinity, made using the local charts U i and V i , i = 1, 2, 3, restricted to the invariant plane z 3 = 0, we have a global picture of the dynamical behavior of system (2) on the sphere S 2 of the infinity, for different values of the parameters a and b, which are shown in Fig. 8.
The integrable systems obtained considering a = 0 or b = 0 are slightly different, but they can be obtained directly from the analysis made above. The phase portraits in the Poincaré sphere of these cases are shown in Fig. 9.
Based on this study at infinity and on the results of Theorem 1 and 2, we can draw a picture of the solutions in the integrable cases (a = 0 and b = 0) contained on the invariant cylinders and on the invariant planes and the ends of these surfaces at infinity, inside the Poincaré ball, as shown in Fig. 10. Observe that the ends of all cylinders at infinity are given by the endpoints of the z-axis, while the ends of all invariant planes z = c are given by the equator of the Poincaré sphere S 2 .
The perturbation of the structures shown in Fig. 10, by making a and b different from zero, and consequently the bifurcation of a periodic orbit from the one contained on the   circle x 2 + y 2 = 2, z = 0, which will be proved in the next section, help us understand the formation of invariant tori in the phase space of Nosé-Hoover oscillator.

The Existence of a Periodic Orbit of System (2) for ab = 0
Based on the study of the integrable cases and on the description of the dynamics of system (2) at infinity, we performed analytical and numerical studies of this system for ab = 0, when it seems to be no more integrable (the non-existence of first integrals of Darboux type in this case was proved in [8]). Using the averaging theory of first order, presented for instance in [2,19], we prove that, for a = 0 and b = 0 sufficiently small, one periodic orbit bifurcates from the circle x 2 + y 2 = 2, z = 0, which is in the family of periodic orbits contained on the invariant plane z = 0 in the integrable case b = 0. More precisely, the following result holds.
Asθ > 0 for ε small enough, we can take θ as the new independent variable in system (11). Doing this and considering the Taylor expansion of order two of the obtained system, we have: Note that differential system (12) is written into the normal form for applying the averaging method (see [2,19]). Using the notation below for system (12): , .

Invariant Tori and Chaotic Behavior of System (2) for ab = 0
The numerical study presented in this section was performed through numerical simulations using the Software Maple™. The equations were solved using a Fehlberg fourth-fifth order Runge-Kutta method with degree four interpolant (known as rk45 method), with step-size equals to 0.01. This method showed to be appropriate for the numerical study of system (2). The files are available to the interested readers, under demand to the authors.

The Formation of Nested Invariant Tori
The periodic orbit γ ε obtained in Theorem 3 bifurcates from the orbit contained in the circle x 2 + y 2 = 2, which belongs to the family of periodic orbits contained in the plane z = 0 in the case b = 0 (see Fig. 10b). This periodic orbit, which persists under small perturbations of the integrable cases, plays an important role in the dynamics of system (2) for a, b different from zero. In fact, by performing a detailed numerical study of system (2) with a and b different from zero, we observe the formation of nested invariant tori around γ ε , which are illustrated in Figs. 11 and 12. Observe that, for ab = 0, only one periodic orbit persists. Also, observe that the solutions which for a = 0 spiral to infinity, upward and downward on the cylinders of radius grater and smaller than 2, respectively, shown in Fig. 10a, connect to each other for ab = 0, forming an invariant torus, as shown in Fig. 11. This happens for all solutions starting near the circle x 2 + y 2 = 2, z = 0, which leads to the creation of nested invariant tori, as shown in Fig. 12.
The proof of the existence of a periodic orbit bifurcating from the circle x 2 +y 2 = 2, z = 0, and the existence of nested invariant tori around it are important from the physical point of view. In fact, as the periodic orbit γ ε is near this circle and the invariant tori are formed around it, in average the position y(t) and the momentum x(t) of the particle modeled by system (2) stay around the circle x 2 + y 2 = 2. Therefore, if we want another average values for these quantities, for instance if we want the solutions staying around the circle Fig. 11 a Periodic orbit of system (2) for a = b = 0.01, and orbits with initial conditions around sea formed through the broken of invariant tori. b Invariant torus, containing inside the periodic orbit γ ε Fig. 12 a Nested invariant tori of system (2) for a = b = 0.01 around the periodic orbit γ ε : taking more initial conditions around γ ε we obtain more invariant tori. b Cross-section showing the existence of nested invariant tori around γ ε x 2 + y 2 = μ 2 , from the proof of Theorem 3 follows that it is enough to consider the following modification of system (2): that is, we introduce a new parameter μ which gives the radius of the circle from which the periodic orbit will bifurcate and, consequently, we will obtain the nested invariant tori around the orbit γ ε,μ , which will lead the solutions to stay close to the circle x 2 + y 2 = μ 2 , z = 0. We observe that the existence of invariant tori for system (1) already appears in the literature as for instance in [13], but not with this approach.

The Formation of Chaotic Behavior
As we vary the parameter values away from the integrable systems (that is, away from a = b = 0), we can observe that some of the nested invariant tori will be broken, giving rise to chaotic sea, as shown in Fig. 13. Other parameter values also lead to the coexistence of   Fig. 15 Graphic of the evolution of Lyapunov exponents for a solution of system (2) with a = 0.6, b = 0.25 and initial condition (x 0 , y 0 , z 0 ) = (0, 0.1, 0). 1.000 iterates were calculated to corroborate the chaotic dynamics of the solutions of system (2) shown in Fig. 13, we calculated the Lyapunov Exponents [21] for a solution with initial conditions (x 0 , y 0 , z 0 ) = (0, 0.1, 0), considering a = 0.6 and b = 0.25, and obtained, after 1.000 iterations, the values: λ 1 = 0.005088, λ 2 = 0.204536, λ 3 = −0.084892, which characterize the chaotic dynamics (see Fig. 15). The occurrence of nested invariant tori coexisting with chaotic motions has already been described in the literature [10,15,16]. It implies the coexistence of a mixing between regions of conservative dynamics with regions of dissipative ones, showing the wealth of such type of differential systems. Observe that system (1) and system (2) with a = 0 have no null divergence, although they present conservative like chaotic dynamics. This subject is treated in the recent published paper [5].

Concluding Remarks and Comments
In this paper, we study globally the dynamics of system (2), which is obtained from the Nosé-Hoover system by rescaling of variables. In order to study and better understand the rich dynamics of Nosé-Hoover system, we consider it as a perturbation of an integrable differential system, obtained considering a = 0 or b = 0 in system (2). We described globally the dynamics of system (2) in the integrable cases, showing the existence of an infinity of concentric invariant cylinders (for a = 0) and an infinity of parallel invariant planes (for b = 0). Using the averaging theory, we proved the existence of a linearly stable periodic orbit which bifurcates from the circle x 2 + y 2 = 2, z = 0, which contains a periodic orbit in the integrable cases. The existence of the bifurcating periodic orbit plays an important role in the existence of nested invariant tori and, consequently, in the breaking of such tori, leading to the formation of chaotic seas surrounding regular islands in the phase space of system (2), for certain parameter values.