Fractal basins of convergence of a seventh-order generalized H\'enon-Heiles potential

This article aims to investigate the points of equilibrium and the associated convergence basins in a seventh-order generalized H\'enon-Heiles potential. Using the well-known Newton-Raphson iterator we numerically locate the position of the points of equilibrium, while we also obtain their linear stability. Furthermore, we demonstrate how the two variable parameters, entering the generalized H\'enon-Heiles potential, affect the convergence dynamics of the system as well as the fractal degree of the basin diagrams. The fractal degree is derived by computing the (boundary) basin entropy as well as the uncertainty dimension.


Introduction
It is well known that every differentiable symmetry of the action of a physical system has a corresponding conservation law. Therefore, by Noether's theorem in every stationary axisymmetric system, the energy and the angular momentum along the symmetry axis are conserved. However, at the end of the XIX century it was shown that in some cases there exists an additional hidden conserved quantity (see e.g. [1,2]), the so-called third integral of motion. This discovery increased the interest of researchers who initiated systematic studies in this topic, among whom Contopoulos stands out by his studies on the existence of the third integral of motion in galactic dynamics [3][4][5][6][7].
An important landmark on the existence of the third integral of motion in axisymmetric potentials is provided by the work of Michael Hénon and Carl Heiles [8], who performed a systematic and complete numerical investigation on this topic, finding that the third integral exists for only a limited range of initial conditions. The potential selected for the study in Ref. [8], can be considered a particular case of the general Hamiltonian found by Contopoulos in [3], H = 1 2 ẋ +ẏ + ω 2 1 x 2 + ω 2 2 y 2 + xy 2 + x 3 .
setting ω 1 = ω 2 = = 1 and = −1/3, and by swapping variables (x, y) → (y, x). The Hamiltonian presented above (1) (and consequently the Hénon-Heiles potential) can be derived as a series expansion up to the third-order of the effective potential for * Corresponding author Email address: evzotos@physics.auth.gr (Euaggelos E. Zotos) stationary axisymmetric systems with reflection symmetry V (r, z) = V (r, −z) (see [9]) Since then, some efforts have been made to generalize the Henon-Heiles potential. Around 1980, Verhulst [10] expanded the potential (2) up to the fourth-order seeking to study resonances 1:1, 1:2, 1:3 and 2:1. Some years ago, a generalized Hénon-Heiles potential was derived by expanding up to the fifth-order the effective potential, aiming to study the equilibrium points and basins of convergence of the new potential [11] and to analyze the dynamical effect on bounded and unbounded orbits of including higherorder terms in the series expansion [12]. More recently, a seventh-order version of the stationary axisymmetric potential was presented [13], finding that when higher-order contributions of the potential are taken into account, the chaoticity of the system is reduced in comparison with the lower-order version of the Hénon-Heiles system. The practical importance of the Henon-Heiles like potentials lies in its applications to the stellar kinematics and velocity ellipsoid in our galaxy, where the observed distribution of star's velocities near the Sun can be explained if a third integral exists [14]. Also, these potentials have been used to investigate quantum manifestations of chaos and level repulsion in classical chaotic Hamiltonians [15], and to calculate the lifetimes and energies for metastable states exploiting the property that the dynamics of this potential changes from quasiperiodic to chaotic for higher energies [16]. In the context of general relativity, these potentials have been used to analyze the emission of gravitational waves and to show the differences among wave emissions from regular and chaotic motion [17], to study the geodesic motion of test particles in vacuum gravitational pp-wave spacetimes [18], or to perform numerical investigations re-lated to the integrability of orbits of test particles moving around a black hole representing the galactic center [19], just to name some examples.
In this paper, we rewrite the general form of the seventhorder potential [13] in terms of two arbitrary parameters α and δ denoting the contributions of the fifth and seventhorder terms, in which the constants are set in such a form that the new potential exhibits an increasing number of fixed points for some values of the free parameters 1 . Aiming to perform a full numerical analysis of the new potential, we shall investigate the existence of equilibrium points using the standard Newton-Raphson iterative scheme. In particular, we will use the so-called basins of convergence [20] in order to explore the optimal initial conditions for which the numerical method is faster and accurate (see e.g. [21][22][23][24]). Moreover, using the probability density function we shall analyze the influence of the free parameters on the convergence of the Newton-Raphson scheme. The fractal degree of the basin diagram will be investigated through the basin entropy and the boundary basin entropy introduced recently by Daza et. al [25][26][27].
The present paper is organized as follows: In section 2, the derivation of the generalized potential along with the new approximate potential is presented. Applying the standard linear stability analysis, in section in 3 the existence and stability of the libration points of the system are calculated as a function of two parameters α and δ related to the contribution of higher-order terms. In section 4, the Newton-Raphson basins of convergence are presented using color code diagrams. Also, we show the biparametric evolution of the basin entropy, the boundary basin entropy and the uncertainty dimension as a function of α and δ. Finally, in section 5 we present the main conclusions of our numerical study.

The model potential
As already pointed out in the introduction section, in a previous paper [13] we derived a generalization of the Hénon-Heiles potential through a Taylor series expansion up to the seventh-order of a generic potential with axial and reflection symmetries. The effective potential is of the form V (r, z) = U (r, z) + L 2 /2r 2 , where r and z denote the radial distance and height of the usual cylindrical coordinates, with V (r, z) = V (r, −z). The seventh-order approximate potential can be written as 1 Note that in Ref. [13] the number of fixed points is always four.
with ξ = r − r 0 , where | * denotes evaluation at (r 0 , 0). It should be pointed out that unlike our previous study, here we redefine the constant factors of the polynomial in order to obtain a large spectrum of fixed points. Also, we introduce two arbitrary parameters α and β, such that setting α = δ = 0 the new potential reduces to the wellknown classical Hénon-Heiles potential. The specific replacements are as follows: Therefore, after applying the previous replacements into Eq. (3), the final potential reads as In the next sections, the main properties and characteristics of the new seventh-order potential are analyzed.

Equilibrium points
The number of points of equilibrium is a function of the values of the parameters α and δ. Our analysis suggests that when α ∈ [0, 10] and δ ∈ [0, 10] we have six cases, depending on the total number of libration points. In Fig. 1 we present the color basins on the (α, δ)-plane which correspond to a different number of points of equilibrium. It is interesting to note, that in all cases the system has always an even number of libration points. Moreover, it is observed that the amount of equilibria becomes mainly affected by the two parameters (α, δ) since the basins do not form vertical or horizontal bands. Fig. 2 shows the equilibrium positions, for eight cases, with values of α and δ, corresponding to all possible combinations of libration points. The coordinates of the libration points are presented as the intersection points of the curves V x = 0 (green lines) and V y = 0 (blue lines). We should note, that in Fig. 1 we have seen that there exist two basins corresponding to 12 points of equilibrium. It turns out that the geometry of the curves V x = 0 and V y = 0, as well as the locations of the equilibrium points, are different in each case. Therefore, we have eight different cases (counting also the classical HH system with α = δ = 0), regarding the total number of libration points.
Once the coordinates of the equilibrium conditions (x 0 , y 0 ) are determined, one can also study their linear stability. The linear stability or instability of a libration point is obtained through the following characteristic equation where V xx , V yy , and V xy denote the second-order partial differentials of the potential V (x, y) with respect to the subindex variable. When the quartic equation (6) has four pure imaginary roots, then the respective point of equilibrium is linearly stable. The existence of four pure imaginary roots is se-cured by the three conditions which must simultaneously be fulfilled.
Our computations indicate the following: • When 4 equilibria exist, only L 1 is linearly stable, while the rest of them are linearly unstable.
• When 6 equilibria exist, only L 1 and L 5 are linearly stable, while the rest of them are linearly unstable.
• When 8 equilibria exist, only L 1 , L 7 , and L 8 are linearly stable, while the rest of them are linearly unstable.
• When 10 equilibria exist, only L 1 and L 5 are linearly stable, while the rest of them are linearly unstable.
• When 12 equilibria exist (the case with the middle blue basin in Fig. 1), only L 1 , L 11 , and L 12 are linearly stable, while the rest of them are linearly unstable.
• When 12 equilibria exist (the case with the upper blue basin in Fig. 1), only L 1 and L 5 are linearly stable, while the rest of them are linearly unstable.
• When 14 equilibria exist, only L 1 , L 11 , and L 12 are linearly stable, while the rest of them are linearly unstable.
The general conclusion is that the point equilibrium located at the origin with x = y = 0, is always linearly stable, regardless of the particular values of the parameters α and δ.

The Newton-Raphson basins of convergence
Knowing the equilibrium positions of a dynamical system is very important. However, in many cases (including our modified HH system) the coordinates of the libration points cannot be derived analytically. Then, the equilibrium solutions can be derived only by employing numerical methods. One of the easiest ways of solving numerically a system of equations (in our case the coupled system V x = V y = 0) is by using the Newton-Raphson (NR) iterative scheme It is a well-known fact, that the outcomes of any numerical method are influenced by the choice of the starting   conditions. In particular, both the speed and the accuracy of any numerical scheme fully depend on the chosen initial conditions. There exist starting conditions for which the iterator diverges, while there are also exist starting conditions leading to one of the roots of the system. The ideal initial conditions (regarding fast convergence and accuracy) form the so-called NR basins of convergence (NR-BoC). This is exactly the importance of identifying the location of the NR-Boc of a dynamical system.
In panels (a)-(h) of Fig. 3 we present the structure of the NR-Boc on the configuration (x, y)-plane, for the eight different cases, classified in terms of the number of equilibrium points. In all cases, the values of the parameters α and δ are the same as those of the panels of Fig. 2. For our computations, the NR scheme was allowed to perform up to 500 iterations, while the desired accuracy, regarding the (x, y) equilibrium positions, was set to 10 −16 .
From the basin diagrams of Fig. 3, it is observed that many structures on the configuration (x, y) plane are very intrincated. Moreover, some of the NR-BoC have a finite domain, while others extend to infinity. Nevertheless, in all cases, there exist well-defined structures containing ideal starting conditions for the numerical scheme. In Fig. 4 we display color maps showing how the required number of iterations N is distributed on the (x, y)-plane. Furthermore, in Fig. 5 we provide the probability distributions. The histograms displayed in Fig. 5 with the probability distributions, may provide additional information about the properties of the modified NR method. For example, the right-hand side of the histograms can be fitted by using the well-known Laplace distribution or double exponential distribution, which is the simplest and most suitable choice [28][29][30].
The probability density function (PDF) for the double exponential distribution reads as where the quantities d > 0 and l are known as the diversity and the location parameter, respectively. Since we are interested only in the probability tails for the histograms, we need only the N ≥ l part of the PDF. We aim to understand how the parameters α and δ influence the convergence properties of the NR scheme. To this end, we defined a 1024 × 1024 grid of (α, δ) values and for each pair, we used the NR scheme for classifying a set of 300 × 300 (x 0 , y 0 ) initial conditions, on the configuration plane and in particular inside the squared region −5 ≤ x, y ≤ +5.
In part (a) of Fig. 6, we present the evolution of the average number of iterations N , needed by the NR method for providing the coordinates of the equilibria with the desired accuracy. Panels (b) and (c) of Fig. 6 depict the distributions of the location parameter (l) and the diversity (d) of the Laplace PDF. Our results strongly indicate that the Laplace PDF is an excellent candidate for fitting the probability histograms, if we take into account that the numerical values of N and l are very close |l − N | ≤ 2). Additionally, from the distribution of the diversity d, shown in part (d), we can conclude that the probability histograms are very well-organized around the average value N , since in most of the cases the numerical value of the diversity is relatively low (d < 5). Finally, in panel (d) of Fig. 6, we show how the differential entropy, defined as h = 1 + ln(2d), evolves as a function of the values (α, δ). It is seen, that both quantities d and h have a very similar parametric evolution. If we take into consideration the combined information from all four panels of Fig. 6 we can argue that the NR method works faster when the system has either 4, 10, 12 or 14 points of equilibrium, while when 6 or 8 libration points exist the convergence of the NR scheme is considerably slower.
Previously, in Fig. 3 we have seen that there are certain regions on the plane (x, y), where using the corresponding starting conditions it is very difficult to know beforehand to which point of equilibrium they are going to converge. These regions are composed of a fractal mixture of final states (equilibria) and they are of course the exact opposite of the basins of convergence. In order to obtain quantitative information about the fractal degree of the BCs on the plane (x, y), we shall compute the basin entropy S b [25,27]. This modern tool indicates the fractal degree of a basin diagram by examining its topological properties. In part (a) of Fig. 7 we show the distribution of the numerical values of S b , as a function of (α, δ). Now we can conclude, without any doubt, that when the system has eight points of equilibrium, we encounter the most fractal NR-BoC, while the fractal degree is considerably lower for a higher number of libration points.
Unfortunately, the transition between smooth and fractal boundaries cannot be determined by the basin entropy S b . The main reason for this drawback is that the basin entropy addresses the uncertainty to link a set of initial conditions to its corresponding final states. Therefore, if we are interested in detecting small variations in the basin boundary we must use another indicator, the boundary basin entropy S bb , which was introduced for the first time in 2016 by Daza et.al. [25]. For obtaining the boundary basin entropy, all we have to do is to divide the total entropy between the number of cells that fall in the boundaries of the convergence basins. This tool gives us the possibility to safely conclude if the basin boundary is fractal or not, by using the so-called "log 2 criterion", with the sufficient condition, if S bb > ln 2, then the boundary is certainly fractal. The distribution of the values of S bb , as a function of (α, δ), is given in panel (b) of Fig. 7. We see that when eight points of equilibrium exist the basin boundaries on the (x, y)-plane are always fractal, while on the other hand when the system has only 4 libration points, the basin boundary entropy exhibits the smaller values when compared to the other cases.
Finally, another standard way to measure the level of fractality of a basin diagram is by computing the fractal dimension [31]. At this point, it is important to emphasize that the results obtained with the basin boundary entropy S bb and the fractal dimension D 0 are related but they do not necessarily have to be the same because the first numerical tool allows us to assess easily that some boundaries are fractal, while the second one provides information about the whole basin since the fractal dimension is an intrinsic property of the system [32,33]. In Fig.  8, we present the dependence of the uncertainty dimension D 0 with the parameters α and δ. As usual, when the fractal dimension equals one, the fractality is zero, while if its value tends to 2 it suggests complete fractality of the respective basin diagram. It is seen, that D 0 displays the highest values when eight points of equilibrium exist, while the lowest values are observed for the cases with 10, 12 and 14 libration points. One should certainly note the large similarity on the parametric evolutionary pattern of D 0 with respect to that of the basin entropy S b . This similarity can be explained by considering that these two computer-based analysis techniques are grounded on boxcounting methodologies.

Discussion
In this work we explored, using numerical techniques, the equilibrium points and the convergence properties of the associated basins of convergence, of a seventh-order generalized Hénon-Heiles potential. The Newton-Raphson root method was used for locating the (x, y) coordinates of the points of equilibrium, while their linear stability was also revealed as a function of both parameters α and δ. Modern color-coded plots were deployed for illustrating the convergence basins on the (x, y) plane. Finally, we managed to determine how the parameters α and δ affect both the accuracy and speed of the NR method, while the fractal degree of the respective basin diagrams was estimated by computing the (boundary) basin entropy and the uncertainty dimension.
The routine of the bivariate NR scheme was coded in FORTRAN 77 (see e.g., [34]). For the taxonomy of the starting points on the plane (x, y) we needed, per grid, roughly about 3 minutes using a Quad-Core i7 4.0 GHz CPU. All the plots of the paper have been developed by using the software Mathematica [35].