Dynamical analysis of bounded and unbounded orbits in a generalized Hénon-Heiles system

Dynamical analysis of bounded and unbounded orbits in a generalized Hénon-Heiles system

F. L. Dubeibe fdubeibe@unillanos.edu.co A. Riaño-Doncel angelica.riano@unillanos.edu.co Euaggelos E. Zotos evzotos@physics.auth.gr Facultad de Ciencias Humanas y de la Educación, Universidad de los Llanos, Villavicencio, Colombia Department of Physics, School of Science, Aristotle University of Thessaloniki, GR-541 24, Thessaloniki, Greece

The Hénon-Heiles potential was first proposed as a simplified version of the gravitational potential experimented by a star in the presence of a galactic center. Currently, this system is considered a paradigm in dynamical systems because despite its simplicity exhibits a very complex dynamical behavior. In the present paper, we perform a series expansion up to the fifth-order of a potential with axial and reflection symmetries, which after some transformations, leads to a generalized Hénon-Heiles potential. Such new system is analyzed qualitatively in both regimes of bounded and unbounded motion via the Poincaré sections method and plotting the exit basins. On the other hand, the quantitative analysis is performed through the Lyapunov exponents and the basin entropy, respectively. We find that in both regimes the chaoticity of the system decreases as long as the test particle energy gets far from the critical energy. Additionally, we may conclude that despite the inclusion of higher order terms in the series expansion, the new system shows wider zones of regularity (islands) than the ones present in the Hénon-Heiles system.

Nonlinear dynamics and chaos – Numerical simulations and chaos – Hamiltonian mechanics – Astrodynamics
journal: Physics Letters A

1 Introduction

The Hénon-Heiles Hamiltonian [1] is considered a representative model for time-independent Hamiltonian systems with two degrees of freedom and thus an essential topic in many textbooks on Nonlinear dynamics, see for example [4, 3, 2]. The main reasons for that are its simple analytical form and at the same time its far from trivial dynamics. Such system was originally formulated to shed lights on the question: does an axisymmetric potential admit a third isolating integral of motion? Nowadays it can be considered one of the most cited works in the field of complex systems, where a huge amount of research has been devoted to discriminate between regular and chaotic motion or to study the escape dynamics of orbits, see e.g. [5, 6, 7, 8, 9, 10, 11, 12]. Although its application was first oriented to the field of galactic dynamics, its applications include semiclassical and quantum mechanics [13, 14, 15].

Some modifications to the original Hénon-Heiles system (henceforth HH-system) have been proposed by different authors, e.g, adding dissipation terms [16], introducing white noise [17, 18], or including forcing terms [19, 20], just to name a few. On the other hand, to our knowledge, the only formal derivation of a generalized Hénon-Heiles potential has been carried out by Verhulst [21], who performed a series expansion up to the fourth-order of a potential with axial and reflection symmetries, finding that his potential differs from the one of Hénon-Heiles by the presence of the quartic polynomial . In addition, it was found that as the area enclosed by the zero velocity curve tends to zero an additional isolating integral can generally be derived [21]. The Verhulst potential has been used, for example, to study the orbital structure near the center of a triaxial galaxy with an analytic core [22], the correlation between the Lyapunov exponents and the size of the chaotic regions in the surface of section [23], or the escape regions in a quartic potential [24].

Following the procedure outlined in Refs. [21] and [25], in the present paper we perform a series expansion of an axisymmetric potential up to the fifth-order, such that the resulting Hamiltonian has as particular cases the Hénon-Heiles and Verhulst systems (abusing terminology, all along the paper we will call the new potential generalized Hénon-Heiles system or simply GHH-system). In order to analyze the gradual transition of the dynamics from the HH-system to the GHH-system, we shall introduce a factor multiplying the quartic and quintic polynomial. Since it is a well-known fact that some types of Hamiltonian systems have a finite energy of escape at which the equipotential surfaces give place to exit channels, here we aim to study the dynamics of the GHH-system not only for energy values below the escape energy (bounded orbits) but also above this threshold value (unbounded orbits).

The paper is organized as follows: In section 2 we calculate the fifth-order series expansion of the potential, and then we write down the Hamiltonian of the GHH-system with their respective equations of motion. Next, we calculate the critical values of energy for bounded and unbounded motions as a function of . In section 3 we study the dynamics of bounded orbits through the Poincaré sections and Lyapunov exponents as a function of the total orbital energy and the parameter . Unbounded orbits are analyzed by means of the exit basins and the basin entropy in section 4. Finally, our paper ends with section 5, where the conclusions are presented.

2 Approximate Potential

The most general Hamiltonian in cylindrical coordinates for a test particle in an axisymmetric potential can be written as


with the component of angular momentum about the -axis.

Let us define the effective potential as , which as a minimum at , where


and corresponds to the radius of a circular orbit on the symmetry plane. By imposing reflection symmetry in , the potential must be an even function in such that all its odd-powered derivatives are odd functions. Expanding the effective potential as a Taylor series around up to the 5th-order, we find


where we have used the fact that an odd function is zero at the origin, we omitted the constant terms, and the remaining constants have been set as follows


With the aim to compare the new potential (2) with the one derived by Hénon and Heiles [1], we replace , , and choose the values , , , such that the new Hamiltonian takes the form111As can be easily noted, setting in Eq. (2) we get the well-known Hénon-Heiles system [25].


The equations of motion derived from Hamilton’s equations read as


where and are the canonical conjugate momenta to and , respectively. Since the total energy is conserved ( constant), the orbital motion is restricted to the region


Depending on the value of the parameter , the dynamical system (6-9) has a given number of real fixed points at which




The introduction of the arbitrary parameter in (2) provides a valuable tool for studying the transition from the HH-system to the GHH-system, for this reason, all along this paper we shall consider values of the parameter in the interval .

0 0.166666666666666667 0.166666666666666667
0.1 0.0905432951155776544 0.101272469408731802
0.5 0.0415033469448885181 0.0475015971035151890
1 0.0267918006221522439 0.0304408354292176817
Table 1: Critical values of energy as a function of for bounded motion and unbounded motion with three escape zones .

On the other hand, in order to classify the orbital motion as bounded or unbounded, we find two critical values of energy, and , such that for energies below the threshold energy there exists a closed zero-velocity surface, while for values larger than the system shows three escape zones (see Table 1). It should be clarified that corresponds to the well-known energy of escape, so for values of larger than , there exists at least one escape channel.

3 Bounded Orbits

In the preceding section, it was possible to determine the value of the energy of escape for the GHH-system and hence to establish that the test particle could exhibit bounded motions for certain values of the total energy , which depend on the values of the parameter . In this section, we study the dynamical behavior of bounded orbits in four different cases and 1, aiming to observe the incidence of the additional term in the HH-potential.

The system of equations (6-9), has been solved using a Runge-Kutta-Fehlberg algorithm RKF8(9) with variable time step. This method lets us numerically solve the system of equations once we know the constants , and the initial conditions and . The initial positions are chosen to fit the condition for confined motion, specifically, we set and in all considered cases. To cover the whole phase space region of allowed motion, we used 15 different initial values of in the interval , where corresponds to the largest value of that allows for a real numerical value of . In each case, is determined by the integral of motion . The existence of the constant of motion indicates that the orbital motion takes place in a three-dimensional effective phase space in which the method of Poincaré sections is an adequate tool to characterize the motion between regular and chaotic.

In Fig. 1, we show the Poincaré sections for the GHH-system using different values of and gradually decreasing the value of . It can be seen from panels , , , and of Fig. 1 that values of very close to the critical energy give rise to chaotic trajectories. On the other hand, for intermediate values of energy the system undergoes a transition from chaos to regularity (see panels , , , and ), while in panels , , , and of Fig. 1 we show that no chaos seems to exist for values of much smaller than the critical energy . Additionally, it can be easily noted from the same figures that the chaos is weaker for larger values of , in other words, the large chaotic zones in the Hénon-Heiles case evolve into regular multi-island orbits endowed in a chaotic sea.

Figure 1: Poincaré surface of sections for (a) and , (b) and , (c) and , (d) and , (e) and , (f) and , (g) and , (h) and , (i) and , (j) and , (k) and , and (l) and . In each panel, the system energy varies according to the relation . The initial conditions have been set as , varying in the interval and is determined by the energy conservation. Here, corresponds to the largest value of that allows for a real numerical value of .
Figure 2: Orbital structure of the plane for (a) and , (b) and , (c) and , (d) and , (e) and , (f) and , (g) and , (h) and , (i) and , (j) and , (k) and , and (l) and . In each panel, the system energy varies according to the relation . The color-code is as follows: regular orbits (green), sticky orbits (blue), chaotic orbits (red). (Color figure online)

The classical Poincaré surfaces of section presented in Fig. 1 give us a quick and rather general overview of the orbital structure in the GHH-system. However, if we require a more detailed view regarding the ordered or chaotic nature of the system then we have to perform a much more thorough scan of the available phase space by classifying large sets of initial conditions of orbits. In Fig. 2 we use color-coded diagrams for revealing the orbital structure of the space. In particular, we define dense uniform grids of initial conditions, regularly distributed inside the limiting curve, and we numerically integrate them for 5000 time units. For obtaining the character of the orbits we use the Smaller Alignment Index (SALI) [26], which is defined as




are the alignments indices, while and , are two deviation vectors which initially are orthonormal and point in two random directions. At every time step, each deviation vector is normalized to 1. Therefore, SALI is a dynamical quantity which inform us if the deviation vectors and have the tendency to obtain the same direction, either by coinciding or by becoming opposite. In the case of chaotic orbits the direction of the deviation vectors has the natural tendency to coincide with that of the most unstable nearby normally hyperbolic invariant manifold, which means that SALI tends to zero. On the contrary, in the case of regular orbits it lies on a torus which directly implies that the two deviation vectors and eventually become tangent to that torus, while generally they converge to entirely different directions. Consequently, the SALI fluctuates around a positive value.

Computationally, the nature of an initial condition is determined according to the final value of SALI at the end of the numerical integration. More precisely, if SALI we have the case of a regular orbit, while if SALI the orbit is chaotic. On the other hand, when we have the case of a sticky orbit222Sticky orbit refers to a special type of orbit which behaves as a regular one for a long-time interval and then it finally exhibits its true chaotic nature.. From Fig. 2, it becomes evident that these color-coded diagrams clearly allow us to identify tiny local stability islands as well as weak chaotic layers, which are hard to be seen in a classical Poincaré surface of section.

A common criterion to measure the chaoticity of a dynamical system is to determine the maximum Lyapunov exponent , which can be computed, for example, using the variational method [27] instead of the well-known two-particle approach, due to the fact that the last one could lead to inconsistent values of (see e.g. [28]). In Fig. 3, we validate the results obtained with the Poincaré sections by calculating the average of an ensemble of trajectories, in terms of the energy for different values of the parameter . In this case, the largest Lyapunov exponent is calculated from the solution of the variational equations of the system, i.e., through the variational method. Because of the numerical nature of , the Lyapunov exponents larger than the threshold value are considered chaotic, while the ones below the threshold are considered regular.

Tracking the evolution of the system, by keeping the initial conditions fixed as and , it can be observed in Fig. 3 that larger values of exhibit a smaller value of , which also decay faster with decreasing energy (it should be noted that we have plotted against , where the energy of the system varies according to the relation ). These results are in agreement with those obtained in Fig. 1, where the phase space is filled with regular islands for small energies.

Figure 3: Energy dependence of in four different cases . In each case, the system energy varies in the interval according to the relation . (Color figure online)

4 Unbounded Orbits

Figure 4: Exit basins of the GHH-Hamiltonian in the configuration space for: (a) and , (b) and , (c) and , (d) and , (e) and , (f) and , (g) and , (h) and , (i) and , (j) and , (k) and , and (l) and . In each panel, the system energy varies according to the relation . The color code is as follows: non-escaping orbits (green); escape through channel 1 (blue); escape through channel 2 (yellow); escape through channel 3 (red). The black circle denotes the scattering region. (Color figure online)
Figure 5: Basin Entropy as a function of the number of sub-matrices in four different cases: (a) , (b) , (c) and (d) . In each panel, the system energy varies according to the relation . (Color figure online)

In this section, we study the trajectories in the GHH-system for energy values larger than the critical energy , in which the zero velocity surfaces are open with three escape channels. To do so, we solve numerically the system of differential equations (6-9) using a RKF-8(9) routine, with 90000 initial positions uniformly distributed along a square grid. Following the usual convention for the escape basins of the HH-system, the initial momenta are determined from the conditions and , with .

All trajectories are classified either into non-escaping or escaping orbits, where the last ones are subdivided according to the exit channel. As noted above, each trajectory can escape to infinity through three different exits, so we use the following convention: exit 1 , exit 2 , and exit 3 . In Fig. 4 we show the exit basin diagrams for the same values of the parameter used in section 3. In this plot, each initial condition is colored according to the escape channel through which it exits, i.e., escape through channel 1 with blue color, escape through channel 2 with yellow color and escape through channel 3 with red color. On the other hand, the green regions denote initial conditions of non-escaping orbits.

Following the same reasoning used in section 3 to trace the evolution of the system, in Fig. 4 we plot the exit basin for the GHH-system using different values of and gradually increasing the value of . The first row (panels , and ) corresponds to the HH-system . In this case, the obtained results are in agreement with the literature [8, 12], since the exit basins become smoother and well-defined as the energy increases. The same behavior is observed for (panels , and ), (panels , and ), and (panels , and ). However, apparently the exit basins are smoother for larger values of . In order to obtain conclusive results, the basin entropy for the invariant sets of the GHH-system will be thoroughly analyzed in the next paragraphs.

The basin entropy was introduced very recently [29] as a new tool to quantitative measure the uncertainty of the basins. Here, the term uncertainty is understood as the difficulty to determine the final state to which a given initial condition will tend to. Unlike all the usual quantities used in nonlinear dynamics (Kolmogorov-Sinai entropy [31, 32], the topological entropy [33], or the expansion entropy [34]), the basin entropy refers to the topology of the basins instead of the evolution of the trajectories itself. For the sake of completeness, let us briefly describe the method for the calculation of the basin entropy.

Let us assume that our dynamical system has attractors (or final states) in a certain region of the phase space, where can be subdivided into a grid composed of square boxes with trajectories per box. Each box contains between 1 and final states, such that we can denote as the probability that inside a box the resulting final state is . Because the trajectories inside a box are independent, the Gibbs entropy of every box can be written as


with the number of final states inside the box .

The entropy of the whole region can be calculated as the sum of entropies of the resulting boxes of the grid, , thus the entropy relative to the total number of boxes (or basin entropy ) is given explicitly by the expression333For a detailed explanation of the method, we refer the interested reader to [29].


Since the final value for basin entropy strongly depends on the total number of boxes (such that for larger values of a more precise value of is obtained), here we use the approach presented in [30] for the calculation of the basin entropy, in which squared-boxes are randomly selected in the phase space region through a Monte Carlo procedure. Following this method, in Fig 5 we have computed the basin entropy as a function of the number of boxes for the exit basins presented in Fig 4. By choosing , it can be observed that in all cases the basin entropy tends to a constant value for .

In panels and of Fig 5 we compute the basin entropy for and 1, respectively. Additionally, in each panel the system energy varies as a function of according to the relation , i.e., we gradually increase the energy values starting very close to the respective critical energy . In analogy with the results of the previous section, it can be observed that as the energy gets closer to the critical energy the basin entropy takes a larger value. Also, it is observed that larger values of significantly reduce the values of the basin entropy, which means that the uncertainty of the basins is larger for the HH-system than for the GHH-system.

5 Conclusions

In the present paper, we propose a new dynamical system that has as particular cases the Hénon-Heiles and Verhulst Hamiltonians. As a key feature, our model has a finite energy of escape which allowed us to study the dynamics of bounded and unbounded orbits with three channels of escape. Via the Poincaré sections method and validated with the Largest Lyapunov exponent, we have shown that bounded orbits are mainly regular in all cases if the total orbital energy is much smaller than the escape energy. In the same vein, unbounded orbits were analyzed through the exit basin structure in configuration space and measuring the basin entropy. Here, we find that the exit basins become smoother and well-defined for larger values of the total orbital energy. Additionally, our numerical investigation suggests that the level of chaos as well as the uncertainty not only depend on the total orbital energy but also on the contribution of higher order terms in the series expansion, i.e, the orbital structure is much more regular in our model than in the one of Hénon-Heiles, indicating that a third integral of motion seems to exist for energy values closer to the escape energy.

Finally, it deserves mentioning that the fractal structure of the basins was also analyzed by means of the boundary basin entropy . Our findings show that in all the considered cases the boundary basin entropy exhibits a similar tendency to the one observed for the basin entropy, such that the sufficient condition for the existence of fractal basin boundaries, , is always satisfied. This result is in agreement with previous studies indicating that these types of fractal structures appear not only in the paradigmatic Hénon-Heiles system but also in a wide variety of open Hamiltonians with two degrees of freedom and three or more escape channels. We expect that our results will be useful in a wide variety complex systems studies, especially those related to the search of a third integral of motion in an axisymmetric potential.


Our thanks go to the anonymous referee for their appropriate suggestions and for carefully reading the manuscript. One of the authors (FLD) gratefully acknowledges the financial support provided by Universidad de los Llanos and COLCIENCIAS, Colombia, under Grant No. 8840.


  • [1] Hénon, M., Heiles, C., Astron. J. 69 (1964) 73.
  • [2] Gutzwiller, M. C., Chaos in classical and quantum mechanics. Vol. 1. Springer Science & Business Media 2013.
  • [3] Hilborn, R. C., Chaos and nonlinear dynamics: an introduction for scientists and engineers. Oxford University Press on Demand 2000.
  • [4] Tabor, M., Chaos and integrability in nonlinear dynamics: an introduction. Wiley 1989.
  • [5] Aguirre, J., Viana, R. L., Sanjuán, M. A. F., Rev. Mod. Phys. 81 (2009) 333.
  • [6] Barrio, R., Blesa, F., Serrano, S., New J. Phys. 11 (2009) 053004.
  • [7] Conte, R., Musette, M., Verhoeven, C., J. Nonlinear Math. Phys. 12 (2005) 212.
  • [8] De Moura, A. P. and Letelier, P.S., Phys. Lett. A 256 (1999) 362.
  • [9] Fordy, A. P., Physica D 52 (1991) 204.
  • [10] Wojciechowski, S., Phys. Lett. A 100 (1984) 277.
  • [11] Zotos, E. E., Nonlinear Dynamics 79 (2015) 1665.
  • [12] Zotos, E. E., Meccanica 52 (2017) 2615.
  • [13] Jaffé, C. and Reinhardt, W. P., The Journal of Chemical Physics 77 (1982) 5191.
  • [14] Feit, M. D., Fleck Jr, J. A., The Journal of chemical physics 80 (1984) 2578.
  • [15] Hamilton, I. P., Light, J. C., The Journal of chemical physics 84 (1986) 306.
  • [16] Seoane, J. M., Sanjuán, M. A. F., Lai, Y. C., Phys. Rev. E 76 (2007) 016208.
  • [17] Seoane, J. M., Sanjuán, M. A. F., Phys. Lett. A 372 (2008) 110.
  • [18] Seoane, J.M., Sanjuán, M.A.F., Int. J. Bifurcat. Chaos 9 (2010) 2783.
  • [19] Blesa, F., Seoane, J. M., Barrio, R. and Sanjuán, M. A.F., International Journal of Bifurcation and Chaos 22 (2012) 1230010.
  • [20] Coccollo, M., Seoane, J.M., Sanjuán, M.A.F., Commun. Nonlinear Sci. Numer. Simul. 18 (2013) 3449.
  • [21] Verhulst, F., Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 290 (1979) 435.
  • [22] De Zeeuw, T., Monthly Notices of the Royal Astronomical Society 215 (1985) 731.
  • [23] Cleary, P. W., Journal of mathematical physics 30 (1989) 689.
  • [24] Barbanis, B., Celestial Mechanics and Dynamical Astronomy 48 (1990) 57.
  • [25] Contopoulos G., Order and Chaos in Dynamical Astronomy, pp. 435. Springer, Berlin 2002.
  • [26] Skokos, Ch., Journal of Physics A: Mathematical and General 34 (2001) 10029.
  • [27] Contopoulos G., Galgani, L. and Giorgilli, A., Phys. Rev. A 18 (1978) 1183.
  • [28] Dubeibe, F. L. and Bermúdez-Almanza, L. D., Int. J. Mod. Phys. C 25 (2014) 1450024.
  • [29] Daza, A., Wagemakers, A., Georgeot, B., Guéry-Odelin, D. and Sanjuán, M. A. F., Scientific reports 6 (2016) 31416.
  • [30] Daza, A., Georgeot, B., Guéry-Odelin, D., Wagemakers, A. and Sanjuán, M. A. F., Phys. Rev. A 95 (2017) 013629.
  • [31] Kolmogorov, A. N., Doklady of Russian Academy of Sciences 119 (1959) 861.
  • [32] Sinai, Y. G., Doklady of Russian Academy of Sciences 124 (1959) 768.
  • [33] Adler, R. L., Konheim, A. G. and McAndrew, M. H., Trans. Amer. Math. Soc. 114 (1965) 309.
  • [34] Hunt, B. R. and Ott, E., Chaos 25 (2015) 97618.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description