This is a preprint of: “Global dynamics and bifurcation of periodic orbits in a modified Nosé- Hoover oscillator”, Jaume Llibre, Marcelo Messias, Alison C. Reinol, J. Dyn. Control Syst., 2020. DOI: [10.1007/s10883-020-09491-5] Global dynamics and bifurcation of periodic orbits in a modified Nosé-Hoover oscillator Jaume Llibre · Marcelo Messias · Alisson C. Reinol Abstract We perform a global dynamical analysis of a modified Nosé-Hoover oscillator, obtained as the perturbation of an integrable differential system. Us- ing 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 av- eraging 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 ex- istence 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. Keywords Nosé-Hoover oscillator · first integral · periodic orbit · averaging theory · invariant tori · chaotic dynamics. Mathematics Subject Classification (2010) 37G15 · 34C25 · 34C37 Jaume Llibre Universitat Autònoma de Barcelona (UAB), Departament de Matemàtiques, 08193, Bel- laterra, Barcelona - CT, Spain. E-mail: jllibre@mat.uab.cat Marcelo Messias Universidade Estadual Paulista (UNESP), Departamento de Matemática e Computação, 19060–900, Presidente Prudente - SP, Brazil. E-mail: marcelo.messias1@unesp.br Alisson C. Reinol Universidade Tecnológica Federal do Paraná, UTFPR, Departamento Acadêmico de Matemática, 86812-460 Apucarana, Paraná, Brazil. E-mail: alissonreinol@utfpr.edu.br 1 Introduction Since its appearance in 1984, the now well-known Nosé-Hoover oscillator, given by the system of ordinary differential equations ẋ = −y − xz, ẏ = x, ż = α (x2 − 1), (1) was widely studied. Obtained from the propositions made by Nosé in the sem- inal 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,12], 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 har- monic oscillator obtained using Nosé’s canonical equations of motion. In this way the time-dependent variables x, y and z represent respectively, the mo- mentum, the position coordinate and the friction coefficient of a particle. As the role of the friction coefficient z is to maintain the average temperature equals 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–14, 16,17,19] and references therein. For details on the physics background for obtaining system (1), see [11] and the nice introduction and motivation section in [12]. More recently, in [14] 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, 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 [17] 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 [17] 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 con- sider 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 ẋ = −y − a xz, ẏ = x, ż = b (x2 − 1), (2) where x, y, z ∈ R3 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 7→ 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 6= 0 is topologically equivalent to system (1). The apparent artificial change of variable considered above to obtain sys- tem (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 anal- ysis of system (2) for a = 0 (or b = 0) and then study what happens for a, b 6= 0 small, that is, we consider perturbations of the integrable cases. 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 is a fundamental dynamical element for the existence of nested invariant tori of system (2) (and, conse- quently, 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 R3 . In sec- tion 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 parame- ter values are varied. Finally, in section 6 we present some concluding remarks and comments. 2 Dynamics of system (2) in the integrable cases Differential system (2) is completely integrable for: a = b = 0; a = 0 and b 6= 0; a 6= 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 equilibrium points, which are centers on the invariant planes z = c, as it is shown in Figure 1. For a = 0 and b 6= 0, system (2) has two independent first integrals, given by H1 (x, y, z) = x2 + y 2 , (3) z y x Fig. 1 Dynamics of system (2) with a = b = 0: the phase space is foliated by invariant planes z = c, c ∈ R, each one having a center at the equilibrium point (0, 0, c). and y  F1 (x, y, z) = 2z − bxy + b(x2 + y 2 − 2) arctan , (4) x with x 6= 0. For b = 0 and a 6= 0, system (2) has the first integral H2 (x, y, z) = z and    2y + axz 2az arctan √  x 4 − a2 z 2  F2 (x, y, z) = (x2 + y 2 + axyz) exp   √ ,  (5) 4 − a2 z 2 with x 6= 0 and z 6= ±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. Theorem 1 Consider system (2) with a = 0 and b > 0. In this case, the following statements hold. (a) The z-axis is invariant under the flow of system (2) and the solution con- tained in it goes from plus to minus infinity. (b) The phase space of system (2) is foliated by the invariant cylinders x2 +y 2 = c, for c > 0, surrounding the z-axis. Moreover, (i) The invariant cylinder x2 + y 2 = 2 is filled by periodic orbits of system (2). (ii) For c > 2, the solutions spiral upwards on the cylinders x2 +y 2 = c. (iii) For 0 < c < 2, the solutions spiral downwards on the cylinders x2 + y 2 = c. If a = 0 and b < 0 the same results hold, reversing appropriately the sense of the solutions in the z-axis and on the invariant cylinders. Proof If a = 0, b > 0 and x = y = 0, then system (2) reduces to ẋ = 0, ẏ = 0, ż = −b. 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 H1 and F1 given by (3) and (4), respectively. Hence the orbits are contained in the intersection {H1 = c} ∩ {F1 = k}, with c, k ∈ R. For c = 2, we have from (4) that the solutions are contained p in the surfaces 2z − bxy = k, k ∈ R. Now, using (3) we have x = ± 2 − y 2 , from which follows that 1 p  z = ± by 2 − y 2 + k , 2 which gives a closed curve for each value of k ∈ R, see Figure 2, which was obtained considering b = 1 and k = 0. Thus, the invariant cylinder x2 + y 2 = 2 is filled by periodic orbits of system (2), as shown in Figure 3 right. This proves statement (b) (i) of Theorem 1 p Fig. 2 Left: Projection of the closed curve z = ± 21 y 2 − y 2 on the yz-plane. Right: Cor- responding periodic orbit of system (2) with a = 0 on the cylinder x2 + y 2 = 2. If x2 + y 2 = c with 0 < c 6= 2, then the orbits are contained in the intersection of the surfaces x2 + y 2 = c and F1 (x, y, z) = k, determined by the curves ! p y 2 2z ± by c − y + b(c − 2) arctan p = k, c − y2 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 x2 + y 2 = c, see Figure 3 left. Statements (b) (ii) and (iii) of Theorem 1 are proved.  Fig. 3 Dynamics of system (2) with a = 0: spiraling behavior of the solutions on the cylinders x2 + y 2 = c, c 6= 2, and the periodic orbits on the cylinder x2 + y 2 = 2. 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). (b) The phase space is foliated by the invariant planes z = c, with c ∈ R. Moreover, in each plane z = c there is only one equilibrium point at (0, 0, c), and we have: (i) If c = 0, then (0, 0, c) is a linear center; (ii) If 0 < c < 2/a, then (0, 0, c) is a stable focus; (iii) If −2/a < c < 0, then (0, 0, c) is an unstable focus; (iv) If c = 2/a, then (0, 0, c) is a stable improper node; if c = −(2/a), then (0, 0, c) is an unstable improper node. (v) If c > 2/a, 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 R3 is foliated by the invariant planes z = c, with c ∈ R. In each invariant plane, system (2) reduces to the following linear differential system ẋ = −y − acx, ẏ = x, (6) where c ∈ R and a > 0. The eigenvalues of the linear part of system (6) at the origin are √ −ac ± a2 c2 − 4 λ= , 2 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. 3 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 one (the Poincaré ball), whose interior is diffeomorphic to R3 and its boundary, the 2–dimensional sphere S2 = {(x, y, z) : x2 + y 2 + z 2 = 1}, plays the role of the infinity of R3 . In this way it is possible to lead the unbounded phase space R3 of system (2) to a compact phase space given by the closed ball of radius one. 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 Ui and Vi , i = 1, 2, 3, with coordinates (z1 , z2 , z3 ) which cover the sphere as a differential manifold. All the points on the invariant sphere (at infinity) in the coordinates of any chart Ui and Vi are obtained simply by making z3 = 0 in the transformed systems. The points in the interior of the Poincaré ball, which is diffeomorphic to R3 , are given in the local charts Ui by z3 > 0 and in the local charts Vi by z3 < 0. See Figure 4 for an illustration of the sphere at infinity S2 and for the orientation of the local charts Ui and Vi . 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 U1 , after making z3 = 0 is ż1 = az1 z2 , ż2 = az22 + b. (7) The z2 –axis is invariant under the flow of system (7). If ab < 0, then system (7) has the equilibrium points r ! b P± = 0, ± − , a p and the straight lines z2 = ± −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 Figure 5 (a). 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 Figure 5 (b). Fig. 4 Local charts Ui , i = 1, 2, 3, used to draw the phase portrait of system (2) on the Poincaré ball. The charts Vi , i = 1, 2, 3, are diametrically opposed to Ui . z2 z2 z1 z1 (a) (b) Fig. 5 Dynamics of system (2) at infinity in the local chart U1 , (a) for a < 0 < b, and (b) for ab > 0. Now we consider the compactification in the chart U2 . After doing z3 = 0 we get the following differential system at infinity ż1 = −az1 z2 , ż2 = bz12 . (8) System (8) has a straight line filled by equilibrium points in the z2 -axis. If b < 0 < a, then the equilibrium points (0, z2 ) are normally hyperbolic, being stable for z2 > 0 and unstable for z2 < 0. Rescaling the system by z1 6= 0, we have that (0, 0) is of saddle type, see Figure 6 (a). For ab > 0 the solutions are contained in the invariant ellipses z12 z2 + 2 = k, a b with k > 0 constant. These ellipses are formed by heteroclinic orbits to the equilibrium points in the z2 -axis, as shown in Figure 6 (b), 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. z2 z2 z1 z1 (a) (b) Fig. 6 Dynamics of system (2) at infinity in the local chart U2 , (a) for b < 0 < a, and (b) for a > 0 and b > 0. In the local chart U3 , after doing z3 = 0, we obtain the following system at infinity ż1 = −z1 (a + b z12 ), ż2 = −bz12 z2 . (9) Note that the z2 -axis is filled by equilibrium points of system (9) and the z1 - axis is invariant under the flow of this system. Moreover, for ab < 0, system (9) has the equilibrium points  r  a P± = ± − , 0 , b while, for ab > 0, there are no equilibrium points outside the z2 -axis. Studying the eigenvalues of the Jacobian matrix at these equilibrium points and rescaling system (9) by z1 6= 0, we obtain the phase portraits shown in Figure 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. z2 z2 z1 z1 (a) (b) Fig. 7 Dynamics of system (2) at infinity in the local chart U3 , (a) for b < 0 < a, and (b) for a > 0 and b > 0. The flow in the local charts Vi , for i = 1, 2, 3, restricted to z3 = 0 is the same as the flow in the respective local charts Ui , reversing the time, because the compactified vector field in Vi is the same as the vector field in Ui 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 Ui and Vi , i = 1, 2, 3, restricted to the invariant plane z3 = 0, we have a global picture of the dynamical behavior of system (2) on the sphere S2 of the infinity, for different values of the parameters a and b, which are shown in Figure 8. (a) (b) Fig. 8 Phase portraits of system (2) at infinity, (a) for b < 0 < a, and (b) for a > 0 and b > 0. The symmetric cases have the same phase portraits, with opposite orientation of the orbits. 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 Figure 9. (a) (b) Fig. 9 Phase portraits of system (2) at infinity, (a) for a = 0 and b > 0, and (b) for b = 0 and a > 0. Based on this study at infinity and on the results of Theorem 1 and 2, we can drawn 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 Figure 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 S2 . (a) (b) Fig. 10 Dynamics of system (2) on the invariant cylinders and on the invariant plane z = 0 inside the Poincaré ball with their ends at infinity for the integrable systems, (a) a = 0 and b > 0, and (b) a > 0 and b = 0. The breaking of structure shown in Figure 10 (a) and (b), obtained making a and b different from zero, and the bifurcation of a periodic orbit from the one contained on the circle x2 + y 2 = 2, z = 0, which will be proved in the next section, help us to understand the formation of invariant tori in the phase space of Nosé-Hoover oscillator. 4 The existence of a periodic orbit of system (2) for ab 6= 0 Based on the study of the integrable systems and on the description of the dynamics at infinity, we perform analytical and numerical studies of system (2) for ab 6= 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]). In order to prove the existence of a periodic orbit in this case, we will use the averaging theory of first order, presented for instance in [2,18]. Using such a method we prove, for a 6= 0 and b 6= 0 sufficiently small, that one periodic orbit bifurcates from the circle x2 + 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, we have following result. Theorem 3 For ab 6= 0, there exists 0 < ε0 ≪ 1 such that for all ε ∈ (0, ε0 ) 6 0 and b̃ 6= 0 small) has a linearly system (1) with a = εã and b = εb̃ (for ã = stable periodic orbit γε , which tends to the circle x2 + y 2 = 2, z = 0 when ε → 0. Moreover, γε is linearly stable if ãb̃ > 0 and unstable if ãb̃ < 0. In this last case γε has a local stable manifold and a local unstable manifold formed each one by two topological cylinders. Proof In order to use the averaging method, we will rescale the parameters a and b, taking a = ε ã and b = ε b̃, with ε > 0 small enough and ãb̃ 6= 0, from which system (2) writes ẋ = −y − ε ã xz, ẏ = x, ż = ε b̃ (x2 − 1). (10) Considering then the change to polar coordinates (r, θ), where x = r cos θ, y = r sin θ, system (10) becomes ṙ = −ε ã r z cos2 θ, θ̇ = 1 + ε ã z cos θ sin θ, ż = ε b̃ (r2 cos2 θ − 1). (11) 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 dr = −εã rz cos2 θ + O(ε2 ), dθ (12) dz = εb̃ (r2 cos2 θ − 1) + O(ε2 ). dθ Note that differential system (12) is written into the normal form for ap- plying the averaging method (see [2,18]). Using the notation below for system (12)     r −ãrz cos2 θ x=y= , t = θ, T = 2π, f (θ, x) = , z b̃(r2 cos2 θ − 1) we have Z 2π     1 −ãrz cos2 θ 1 −ãrz f0 (x) = dθ = . 2π 0 b̃(r2 cos2 θ − 1) 2 b̃(r2 − 2) Therefore √ f0 (x) = 0 has a unique solution with r > 0, namely p = (r, z) = ( 2, 0), which satisfies det(Dx f0 (p)) = ãb̃ 6= 0. By the averaging method, this solution correspond to a periodic orbit φ(θ, ε) = (r(θ, ε), z(θ, ε)) √ of system (12), for ε > 0 sufficiently small, such that (r(0, ε), z(0, ε)) → ( 2, 0) when ε → 0. Going back through the changes of coordinates we have x2 + y 2 = 2, from which the statement of the theorem follows. Furthermore, p the eigenvalues √ of the Jacobian matrix (Dx f0 )( 2, 0) are given by λ1,2 = ± −ãb̃, from which follows that the periodic solution is linearly stable if ãb̃ > 0, and unstable of saddle type if ãb̃ < 0. Hence Theorem 3 is proved.  5 Invariant tori and chaotic behavior of system (2) for ab 6= 0 5.1 The formation of nested invariant tori The periodic orbit γε obtained in Theorem 3 bifurcates from the orbit con- tained in the circle x2 + y 2 = 2, which belongs to the family of periodic orbits contained in the plane z = 0 in the case b = 0 (see Figure 10 (b)). This pe- riodic 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 γε . The numerical study presented here was performed through numer- ical 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). In Figures 11 and 12 we illustrate the formation of invariant tori around the bifurcated periodic orbit. First, observe that for ab 6= 0, only one periodic orbit persist. Also, observe that the solutions which for a = 0 spiral to infinity, upward and downward on the cylinders of radius grater and smaller than two, respectively, shown in Figure 10 (a), connect to each other for ab 6= 0, forming an invariant torus, as shown in Figure 11. This happens for all solutions starting near the circle x2 + y 2 = 2, z = 0, which leads to the creation of nested invariant tori, as shown in Figure 12. (a) (b) Fig. 11 (a) In red the periodic orbit γε of system (2) for a = b = 0.01, and in blue orbits with initial conditions around the periodic orbit, showing the formation of an invariant torus. (b) Invariant torus, containing inside the periodic orbit γε . The proof of the existence of a periodic orbit bifurcating from the circle x2 + y 2 = 2, z = 0, and the existence of nested invariant tori around it are (a) (b) 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 γε . important from the physical point of view. In fact, as the periodic orbit γε are 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 x2 + y 2 = 2. On the other hand, if we want another average values for these quantities, for instance if we want the solutions stay around the circle x2 + y 2 = µ2 , from the proof of Theorem 3 follows that it is enough to consider the following modification of system (2): ẋ = −y − a xz, ẏ = x, ż = b (x2 − µ2 /2), (13) 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 x2 + y 2 = µ2 , z = 0. The existence of invariant tori for system (1) already appear in the litera- ture as for instance in [12], but not with this approach. 5.2 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 broken, giving rise to chaotic sea, as shown in Figure 13. Other parameter values also lead to the coexistence of nested invariant tori with chaotic behav- ior, as shown in Figures 13 and 15. It seems from the numerical simulations that the chaotic sea is created starting with the destruction of the more ex- ternal invariant tori in the nested invariant tori structure around the periodic orbit γε . In order to corroborate the chaotic dynamics of the solutions of sys- tem (2) shown in Figure 13, we calculated the Lyapunov Exponents [20] for a solution with initial conditions (x0 , y0 , z0 ) = (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 (se also Figure 14). (a) (b) Fig. 13 (a) Chaotic behavior of system (2) for a = 0.6 and b = 0.25: orbit in a chaotic sea surrounding the nested invariant tori. (b) The 2π-Poincar map shows that there is a chaotic sea near the most external torus. 1 0.5 0 −0.5 0 200 400 600 800 1000 Fig. 14 Graphic of the evolution of Lyapunov exponents for a solution of system (2) with a = 0.6, b = 0.25 and initial condition (x0 , y0 , z0 ) = (0, 0.1, 0). 1.000 iterates were calculated. The occurrence of nested invariant tori coexisting with chaotic motions have already been described in the literature [10,14,15]. It implies the coexis- tence of a mixing between regions of conservative dynamics with regions of dis- sipative ones, showing the wealth of such type of differential systems. Observe Fig. 15 In blue is shown one of the nested invariant tori of system (2) for a = 1 and b = 0.25; in red an orbit in the chaotic sea formed through the broken of invariant tori. The 2π-Poincaré map shows a cross section of the invariant tori and the chaotic sea surrounding regular regions. that system (1) and system (2) with a 6= 0 have no null divergence, although they present conservative like chaotic dynamics. This subject is treated in the recent published paper [5]. 6 Concluding remarks and comments In this paper we study globally the dynamics of system (2), which is obtained from the Nosé-Hoover system by a rescaling of variable. 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 x2 + y 2 = 2, z = 0, which contains a periodic orbit in the integrable cases. The existence of the bifurcating periodic orbit plays a important role in the existence of nested invariant tori and, consequently, in the breaking of such tori, leading to the formation of chaotic sea surrounding regular islands. Acknowledgements The first author was partially supported by a MINECO-FEDER grant MTM2016-77278-P, a MINECO grant MTM2013-40998-P, and an AGAUR grant num- ber 2014SGR-568. The second author was partially supported by CNPq-Brazil under the grant 311355/2018-8. The second and third authors were supported by São Paulo Research Foundation (FAPESP) grants numbers 2013/24541-0 and 2018/23190-3. References 1. Cima A, Llibre J. Bounded polynomial vector fields. Trans Amer Math Soc. 1990; 318:557–579. 2. Guckenheimer J, Holmes P. Nonlinear oscillations, dynamical systems and bifurcations of vectors fields. Springer-Verlag: New York; 1983. 3. Holian BL, Voter AF. Thermostatted molecular dynamics: How to avoid the Toda demon hidden in Nose-Hoover dynamics. Phys Rev E. 1995; 52:2338–2347. 4. Hoover WG. Canonical dynamics: Equilibrium phase-space distributions. Phys Rev A. 1985; 31:1695–1697. 5. Jafari S, Sprott JC, Dehghan S. Categories of Conservative Flows. Int J Bifurcat Chaos. 2019; 29:1950021 (16 pages). 6. Llibre J, Messias M. Global dynamics of the Rikitake system. Physica D. 2009; 238:241– 252. 7. Llibre J, Messias M, da Silva PR. Global dynamics in the Poincar ball of the Chen system having invariant algebraic surfaces. Int J Bifurcat Chaos. 2012; 22:1250154 (17 pages). 8. Mahdi A, Valls C. Integrability of the Nos–Hoover equation. J Geom Phys. 2011; 61:1348– 1352. 9. Messias M. Dynamics at infinity and the existence of singularly degenerate heteroclinic cycles in the Lorenz system. J Phys A: Math Theor. 2009; 42:115101 (18 pages). 10. Messias M, Reinol AC. On the existence of periodic orbits and KAM tori in the Sprott A system: a special case of the Nos-Hoover oscillator. Nonlinear Dynam. 2018; 92:1287–1297. 11. Nosé S. A unified formulation of the constant temperature molecular-dynamics methods. J Chem Phys. 1984; 81:511–519. 12. Posch HA, Hoover WG, Vesely FJ. Canonical dynamics of the Nosé oscillator: Stability, order, and chaos. Phys Rev A. 1986; 33:4253–4265. 13. Rech PC. Quasiperiodicity and Chaos in a Generalized Nosé-Hoover Oscillator. Int J Bifurcat Chaos. 2016; 26:16501701 (7 pages). 14. Sprott JC, Hoover WG, Hoover CG. Heat conduction, and the lack thereof, in time- reversible dynamical systems: Generalized Nosé-Hoover oscillators with a temperature gradient. Phys Rev E. 2014; 89:042914 (6 pages). 15. Sprott JC. A dynamical system with a strange attractor and invariant tori. Phys Lett A. 2014; 378:1361–1363. 16. Sprott JC. Strange attractors with various equilibrium types. Eur Phys J Special Topics. 2015; 224:1409–1419. 17. Swinnerton-Dyer P, Wagenknecht T. Some third-order ordinary differential equations. Bull London Math Soc. 2008; 40:725–748. 18. Vehrulst F. Nonlinear differential equations and dynamical systems. Universitext, Springer-Verlag: Berlin Heidelberg; 1996. 19. Wang L, Yang X-S. The invariant tori of knot type and the interlinked invariant tori in the Nosé-Hoover oscillator. Eur Phys J B. 2015; 88:78 (5 pages). 20. Wolf A, Swift JB, Swinney HL, Vastano JA. Determining Lyapunov exponents from a time series. Physica D 1985; 16:285–317.