On the dynamics of a seventh-order generalized H\'enon-Heiles potential

This paper deals with the derivation and analysis of a seventh-order generalization of the H\'enon-Heiles potential. The new potential has axial and reflection symmetries, and finite escape energy with three channels of escape. Based on SALI indicator and exits basins, the dynamic behavior of the seventh-order system is investigated qualitatively in cases of bounded and unbounded movement. Moreover, a quantitative analysis is carried out through the percentage of chaotic orbits and the basin entropy, respectively. After classifying large sets of initial conditions of orbits for several values of the energy constant in both regimes, we observe that when the energy moves away from the critical value, the chaoticity of the system decreases and the basin structure becomes simpler with sharper and well defined bounds. Our results suggest that when the seventh-order contributions of the potential are taken into account, the system becomes less ergodic in comparison with the classical version of the H\'enon-Heiles system.


Introduction
In 1963 Michael Hénon and Carl Heiles, numerically investigated the existence of the third integral of motion in axisymmetric potentials [1], motivated by the apparent existence of such quantity in a large number of galactic orbits. This system is usually called the Hénon-Heiles Hamiltonian and, to our knowledge, it was the first one for which a third integral of motion was constructed [2]. Although its first application was the explanation of the non-equality on the axes of the velocity ellipsoid for stars in the solar neighborhood, its current uses range from the calculation of energies and lifetimes for metastable states within the quantum mechanics formalism [3], to the analysis of gravitational waves, by means of general relativity [4].
Exploiting the properties of Hamiltonian systems with two degrees of freedom (see e.g. [5]), over the last few decades many real problems coming from Physics and Chemistry have been formulated using modifications to the Hénon-Heiles potential. For example, the chaotic ionization mechanism in chemical reactions is modeled with a rotating potential of HénonHeiles type [6], laser-driven reactions with several open channels are modeled using a driven Hénon-Heiles system [7], the description of motion of test particles in presence of vacuum gravitational pp-wave spacetimes was studied using a modified Hénon-Heiles potential [8], or the modeling of black holes with external halos is realized through the superposition of quadrupolar and octupolar terms, such that in the Newtonian limit, the gravitational interaction reduces to an analog of the Hénon-Heiles potential [9], just to name a few.
On the other hand, many efforts have been made to generalize the Hénon-Heiles potential, keeping its main characteristic intact, i.e., being a time-independent two-dimensional Hamiltonian system with finite escape energy. The very first generalization was performed by Verhulst [10], who extended up to the fourth-order this potential, to produce qualitative and quantitative results for the study of the resonance cases 1:2, 1:1, 2:1 and 1:3. More recently, the Hénon-Heiles potential was extended up to the fifth-order, where the authors studied the basins of convergence of equilibrium points of the new system [11] and carried out a dynamical analysis of unbounded and bounded movement [12].
Our goal in this paper is twofold. First, we shall perform a Taylor series expansion up to the seventh-order of a general form of a potential with reflection and axial symmetries, such that the derived Hamiltonian reduces to the lower order cases: Hénon-Heiles, Verhulst and the fifth-order system. Second, by using some state-of-the-art numerical techniques like SALI indicators and exits basins, the dynamic behavior of the new system is investigated in the cases of open and closed zero-velocity curves. The quantitative analysis is carried out through the percentage of chaotic orbits and the basin entropy, respectively.
The present paper is organized as follows. In section 2, the derivation of the generalized potential along with the Hamiltonian and equations of motion are presented. In section 3, the existence and stability of the fixed points are investigated by means of the standard linear stability analysis. Next, we describe the methodology and techniques used to obtain the numerical results for bounded and unbounded orbits. Also, we show how the configuration (x, y) and the (x, E) planes evolve, as a function of a perturbation parameter. Finally, in section 4 we present the main conclusions and implications for the complex systems field, derived from this work.

Derivation of the Potential
Let us consider the case of a test particle in the presence of an axially symmetric potential. Given the symmetry properties of the system, the usual and appropriate coordinate system of choice is cylindrical (r, θ, z). Moreover, the symmetry with respect to the rotational degree of freedom gives place to a cyclic coordinate and therefore, the Hamiltonian can be written as with Φ eff (r, z) = Φ(r, z) + L 2 z /2r 2 , the effective potential and L z , the angular momentum along the z−axis.
Following Contopolulos [13], it is possible to find an approximate analytical form of the axisymmetric potential by performing a Taylor series expansion of Φ eff (r, z) around (r 0 , 0) up to the 7th-order. Under the following considerations: (i) The effective potential Φ eff (r, z) has a minimum at (r, z) = (r 0 , 0), where ∂Φ/∂r = L 2 z /r 3 , and (ii) By imposing reflection symmetry about the z = 0 plane, the potential must contain only even powers on z and hence the odd- where ξ = r − r 0 , we omitted constant terms, and the terms Let us enumerate the main properties of the potential presented above: 1. It has finite escape energy, such that after some energy threshold the zero-velocity curves (henceforth ZVC) exhibit three channels of escape. 2. By setting δ = 0 in Eq. (3), it reduces to the well known classical Hénon-Heiles potential. 3. With an appropriate setting of the constant parameters, it reduces to the fourth-order [10] and fifth-order [12] generalizations of the Hénon-Helies potential.
On the other hand, given the dependence of the Hamiltonian with the velocities, the equations of motion can be written in compact form as  Taking into account that the system is scleronomous and there are no velocity-dependent terms in the potential energy, we can write H = E, in such a way that the orbits must be restricted to the region satisfying The escape energy E esc , is defined as the smallest energy value for which the ZVC open and the test particle is allowed to escape to infinity. In the present study we are interested in bounded and unbounded motion, however, aiming to compare with the Hénon-Heiles system in which energies above E > E esc allows the test particle to escape to infinity through three different channels, in Table 1 we present the critical values of energy for different values of δ satisfying both requirements, i.e., E min refers to the largest energy value for which the ZVC are closed, while E max defines the smallest energy value for which the ZVC has three escape channels.

Results
We start by considering the fixed points of the system and their stability. The system always exhibits four equilibrium points, regardless of the value of δ, such equilibriums shall be identified as L 1 , L 2 , L 3 , and L 4 . In Fig. 1 we present the open ZVC of the system, along with the positions of the equilibrium points for δ = 0 and δ = 1, respectively. At the same diagram, we also depict the positions of the Lyapunov orbits [14,15], which control the escape of the particles. When δ varies in the interval [0, 1] the fixed point L 1 keeps at the origin (0, 0), while L 2 stays static at (0, 1). Only the equilibrium points L 3 and L 4 change position. In Fig. 2 we show the parametric evolution of the coordinates x(L 3 ) and y(L 3 ), as a function of δ. The evolution suggests that for increasing values of δ the equilibrium points L 3 and L 4 move toward the origin, following a nonlinear path. The positions of x(L 4 ) and y(L 4 ) are not presented because they are symmetric with respect to the y−axis, i.e., x(L 4 ) = −x(L 3 ) and y(L 4 ) = y(L 3 ). Concerning the stability, it is found that all libration points are always linearly unstable, when δ ∈ [0, 1].
In what follows, we analyze the dynamics of the seventh order generalized Hénon-Heiles potential. In general terms, we shall classify the bounded and unbounded orbits by using color coded diagrams to distinguish between the different types of the orbits, i.e. for closed ZVC, the classification corresponds to regular, chaotic or sticky orbits, while for open ZVC, the dynamic behavior is richer adding three additional final states, depending on the exit channels: exit via channel 1, exit via channel 2 and exit via channel 3. Once the classification is made, we will use quantitative indicators to inspect the evolution of chaoticity in comparison with the classical Hénon-Heiles system. Let us start with the case in which energy demands all orbits to be bounded.

Bounded orbits
The dynamics of the system in the case of closed ZVC is studied in the present subsection. Here, aiming to reveal the differences in the dynamical behavior of the seventh-order generalization and the Hénon-Heiles system, we shall use the chaos indicator SALI [16]. In particular, we use different values of δ, allowing us to gradually control the contribution of higher-order terms in the system. Moreover, we decrease the energy with respect to its respective critical value to see the influence on the global dynamics of the system.
For a full overview of the dynamics, we define a dense grid of 1024 × 1024 initial conditions, uniformly distributed and closely spaced on the configuration plane (x, y). The initial conditions are set in such a way that they satisfy the energy restrictions, i.e., for each initial coordinate position (x 0 , y 0 ) withẏ 0 = 0 the initial condition forẋ 0 = 0 is always defined by the relation H = E. All our numerical calculations were performed with a double-precision adaptive Bulirsch-Stoer algorithm, implemented in FORTRAN 77 [17], where it was found that by decreasing the magnitude of the time step to 10 −2 , the numerical error has been reduced to an approximate order of 10 −13 , or smaller. It is important to note that we use SALI for distinguishing between chaotic and ordered motion instead of any other indicator (e.g. GALI, LCN, FLI, etc.) because it is not only fast and easy to compute but also efficient and very accurate. According to the literature, orbits are classified as regular if the value of SALI for a maximum integration time (10 4 in our case) is larger than 10 −4 , or chaotic if the resulting value is smaller than 10 −8 , however, if the final value for SALI belongs to the interval 10 −8 ≤ SALI ≤ 10 −4 the orbit is classified as sticky. In this context, sticky is referred to an orbit that needs a larger integration time in order to reveal its chaotic nature 1 .
In Fig. 3 we classify each orbit according to SALI, for different values of δ and energy E. Each row in Fig. 3 has a different value of δ, from up to down: 0, 0.3, 0.6, 1. Also, each column corresponds to a different value of energy varying according to the expression E = E min (1 − n/100), such that we reduce the total energy in n% with respect to the critical energy, in particular, n takes (from left to right) the values: 1, 5, 10, 20. The color code used for the final value of SALI is presented in the lower part of Fig. 3, where the rainbow color palette ranges from dark blue (SALI = 10 −15 ) to dark red (SALI = √ 2).
When reading the figure from left to right it is found that as the energy becomes smaller than the critical energy E min , the area occupied by the regular regions increases in comparison with the area of the chaotic orbits. This behavior has been studied extensively in the theoretical and numerical regimes for the Hénon-Heiles system (see e.g. [18]), but there seems to be also present for δ ∈ (0, 1]. Furthermore, when reading the figure from up to down, it is observed that increasing the values of the δ parameter reinforces the tendency of the system to increase the area occupied by regular orbits, or in other words, from a qualitative point of view, the amount of chaotic orbits seems to decrease as δ increases. A quantitative measure of the percentage of chaotic orbits is presented in Fig. 4. The proportion of chaotic orbits (relative area covered by the blue pixels) is calculated as a function of the parameter δ, using four energy levels, E 1 , E 2 , E 3 , and E 4 , which correspond to n = 1, n = 5, n = 10, and n = 20, respectively. In the four cases we observe that when δ tends to 1 (additional terms equally weighted), the percentage of chaotic orbits is reduced in approximately 40% of its initial value for δ = 0. Additionally, it can be easily noted that the percentage of regular orbits approximately reach the hundred percent for δ = 1 and E = E 4 . These results suggest that when the seventh order contributions of the potential are taken into account, the system becomes less ergodic in comparison with the classical version of the Hénon-Heiles system and in some cases practically regular.

Unbounded orbits
The dynamic analysis of the system in the case of open ZVC with escape channels is studied in the present subsection. In this case, the exit basins are plotted using energy values beyond the threshold limit E max which directly depends on the δ parameter, as shown in Table 1. In general, exit basins are defined as the set of initial conditions that escape through a particular exit, hence, the complexity of the basins may depend on the number of exits on the ZVC. However, as in the present study the system always exhibits three channels of escape, the structure of the exit basins shall be related uniquely to the effect of the δ parameter and the total energy of the system E.
The initial conditions for each orbit are set by defining a position coordinate r = (x, y), while the initial velocities are defined by r ·ṙ = 0 and E = H, with r ×ṙ > 0. For the classification of the orbits we use the standard convention, i.e., if the test particle escapes through exit i, its initial conditions will belong to the exit i-th basin. The exits present in our system are denoted as follows: a particle escaping to infinity through y → ∞ belongs to the channel 1 (green), if (x → −∞, y → −∞) belongs to the channel 2 (red), or if (x → ∞, y → −∞) belongs to the channel 3 (red). It is a well-known fact that for energy values larger than the escape energy a portion of orbits does not escape from the system, and sometimes forming islands around the stable periodic orbits [19], this phenomenon occurs mainly for energy levels, just above the critical value of the energy of escape.
In Fig. 5 we present the exit basins obtained for four different values of δ and energy. Each row has a different value of δ, from up to down: 0, 0.3, 0.6, 1, while, each column corresponds to a different value of energy varying according to the expression E = E max (1+n/100), such that n corresponds to the percantage increase in energy with reference to the critical value E max and takes the values (from left to right): 1, 5, 10, and 20. The first row corresponds to the Hénon-Heiles system, where it can be noted that for larger percentages of deviation in respect of E max , the basin boundaries become smoother and also well-defined, in fact, for n ≥ 50 the amount of orbits escaping through the channels is so high that the bounded orbits become imperceptible to the naked eye. The same behavior is perceived for δ ∈ (0, 1] (second, third and fourth row), nevertheless, unlike the previous case, the lower exit channels show a symmetric wider bottleneck thus allowing most orbits to go to infinity through these two channels.
The global dynamical behavior of the system can be analyzed by means of the orbit classification in the (x, E) plane. In Fig. 6 we present the final states of the system using a continuous spectrum of values of energy, for six values of δ = 0, 0.2, 0.4, 0.6, 0.8, and 1. The initial conditions on xposition and energy are selected from the allowed regions of motion, while y =ẋ = 0 andẏ is determined by the energy conservation E = H. This figure shows that for larger values of E the phase space is mainly occupied by escaping orbits. For energy values close to E min , practically all the final states are possible but the basin boundaries are poorly defined and the extent of these regions increases with the increase in the parameter δ. Finally, when E → 0 the orbits of the system are predominantly regular with a chaotic layer surrounding it, however, this region becomes less uniform as δ → 1.
On the other hand, a quantitative measure of the uncertainty of the exit basins is provided by the basin entropy [20]. Such an indicator gives us information about the topology of the basins rather than of the evolution of the orbits, or in other words, a measure of the difficulty to predict in advance which will be the final state of an orbit. In a generic case, for a dynamical system with N A final states, whose phase space was subdivided into a numerical grid of N square boxes, the basin entropy can be calculated as the average of the Gibbs entropy of every box, i.e., where P i, j represents the probability to get the j-th final state inside the i-th box 2 .
In Fig. 7 we plot the basin entropy as a function of the δ parameter using four levels of energy, E 1 , E 2 , E 3 , and E 4 , which correspond to n = 2, n = 50, n = 100, and n = 200, respectively. In all cases, we observe that as the parameter δ increases, the basin entropy S b reduces, such decrease is more marked for small δ values and becomes less pronounced (with a slope close to zero) for high δ values. Also, it is important to note, that there is no strict relationship for fixed values of δ and decreasing energies, for example it is observed that when δ = 0 (Hénon-Heiles system) the basin entropy is larger for levels E 4 than for levels E 3 , or even for a full contribution of the seventh order terms δ = 1 the basin entropy is much smaller for levels E 1 than for levels E 2 .

Conclusions
In this paper, by using some state-of-the-art numerical techniques, we have performed a numerical investigation of the dynamics of bounded and unbounded orbits with three channels of escape for a seventh-order generalization of the Hénon-Heiles potential. Our dynamical model can be used, for example, to describe the motion of stars in the vicinity of the galactic center, since it comes from the Taylor series expansion of an axially symmetric galactic potential, with reflection symmetry about the plane z = 0.
A systematic classification of orbits is carried out, so that the different types of orbits can be arranged into groups, according to the structure of the zero-velocity curves. For closed contours the orbits are classified as regular, chaotic or sticky, while for open contours with three channels of escape, orbits are assorted as regular, chaotic or sticky, escaping through exit 1, escaping through exit 2 or escaping through exit 3. The method used to determine the chaoticity or regularity of the orbits is the standard SALI method, which provides a reliable and fast indicator of chaos in Hamiltonian systems.
To obtain a graphical representation of the dynamics, we define a dense grid of initial conditions, uniformly distributed and closely spaced on the configuration plane (x, y). The qualitative information supplied by the basins suggests that, in accordance with the well-known results for the classical Hénon-Heiles Hamiltonian, for energy values larger or smaller than the critical value, the chaoticity of the system decreases and the structure becomes simpler with sharper and well-defined bounds. This fact is confirmed when plot-ting the basins on the (x, E) plane, thus using a continuous spectrum of energy levels.
Concerning the effect of the seventh-order terms in the orbital dynamics of the system, we studied the quantitative evolution of the percentage of chaotic orbits and the basin entropy for the exit basins, as a function of a parameter δ. Our results are conclusive, pointing out that the percentages of chaotic orbits are significantly reduced for the new potential in comparison with the results for the classical Hénon-Heiles, even, for very small energy values, the new system approximately reaches the hundred percent or regular orbits. Similarly, the highly fractal basin boundaries of the classical Hénon-Heiles become smoother and well-defined when a full contribution of the seventh order terms is present. We hope that our results will provide a reference for future research on the field of complex systems, mainly in those studies in which a third integral of motion exists in axisymmetric potentials.