Abstract
This article aims to investigate the points of equilibrium and the associated convergence basins in a seventh-order generalized Hénon–Heiles potential. Using the well-known Newton–Raphson iterator, we numerically locate the positions of the points of equilibrium, while we also obtain their linear stability. Furthermore, we demonstrate how the two variable parameters, entering the generalized Hénon–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.
1. 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 19th century, it was shown that in some cases there existed 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–7].
An important landmark on the existence of the third integral of motion in axisymmetric potentials is provided by the work of Hénon and Heiles [8], who performed a systematic and complete numerical investigation on this topic and found that the third integral exists for only a limited range of initial conditions. The potential selected for the study in [8] can be considered a particular case of the general Hamiltonian found by Contopoulos in [3]:
Set and , and swap variables .
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 stationary axisymmetric systems with reflection symmetry (see [9]).
Since then, some efforts have been made to generalize the Hénon–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 the effective potential up to the fifth order, 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 higher-order terms in the series expansion [12]. More recently, a seventh-order version of the stationary axisymmetric potential was presented [13], and it was found 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 Hénon–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 change 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 related 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 seventh-order potential [13] in terms of two arbitrary parameters and denoting the contributions of the fifth-order and seventh-order 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 (note that in [13] the number of fixed points is always four). 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] to explore the optimal initial conditions for which the numerical method is faster and accurate (see, e.g., [21–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–27].
The remainder of this 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 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.
2. The Model Potential
As already pointed out in the Introduction, 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 , where and denote the radial distance and height of the usual cylindrical coordinates, with . The seventh-order approximate potential can be written aswithwhere denotes evaluation at .
It should be pointed out that, unlike our previous study, in the present paper, we redefine the constant factors of the polynomial with two main objectives: first, we introduce two arbitrary parameters and , such that, setting , the new potential reduces to the well-known classical Hénon–Heiles potential, and, second, we aim to obtain a large spectrum of fixed points. This last process was carried out numerically by varying each coefficient and observing the total number of fixed points until obtaining the combination that exhibits the richest dynamics for the system. The specific replacements are as follows: , , , , , , , , and .
Therefore, after applying the previous replacements into equation (3), the final potential reads
In the next sections, the main properties and characteristics of the new seventh-order potential are analyzed.
3. Equilibrium Points
The number of points of equilibrium is a function of the values of the parameters and . Our analysis suggests that when and , we have six cases, depending on the total number of libration points. In Figure 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 number of equilibria becomes mainly affected by the two parameters since the basins do not form vertical or horizontal bands.
Figure 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 (green lines) and (blue lines). We should note that in Figure 1 we have seen that there exist two basins corresponding to 12 points of equilibrium. It turns out that the geometry of the curves and , as well as the locations of the equilibrium points, is different in each case. Therefore, we have eight different cases (counting also the classical HH system with ), regarding the total number of libration points.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
Once the coordinates of the equilibrium conditions 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 , , and denote the second-order partial differentials of the potential with respect to the subindex variable.
When the quartic equation (6) has four pure imaginary roots, the respective point of equilibrium is linearly stable. The existence of four pure imaginary roots is secured by the three following conditions:which must simultaneously be fulfilled.
Our computations indicate the following:(i)When 4 equilibria exist, only is linearly stable, while the rest of them are linearly unstable(ii)When 6 equilibria exist, only and are linearly stable, while the rest of them are linearly unstable(iii)When 8 equilibria exist, only , , and are linearly stable, while the rest of them are linearly unstable(iv)When 10 equilibria exist, only and are linearly stable, while the rest of them are linearly unstable(v)When 12 equilibria exist (the case with the middle blue basin in Figure 1), only , , and are linearly stable, while the rest of them are linearly unstable(vi)When 12 equilibria exist (the case with the upper blue basin in Figure 1), only and are linearly stable, while the rest of them are linearly unstable(vii)When 14 equilibria exist, only , , and are linearly stable, while the rest of them are linearly unstable
The general conclusion is that the point equilibrium located at the origin with is always linearly stable, regardless of the particular values of the parameters and .
4. 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 ) 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 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 Figures 3(a)–3(h), we present the structure of the NR-BoC on the configuration -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 Figure 2. For our computations, the NR scheme was allowed to perform up to 500 iterations, while the desired accuracy, regarding the equilibrium positions, was set to .
From the basin diagrams of Figure 3, it is observed that many structures on the configuration -plane are very complicated. 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 Figure 4, we display color maps showing how the required number of iterations is distributed on the -plane. Furthermore, in Figure 5, we provide the probability distributions.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
The histograms displayed in Figure 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–30].
The probability density function (PDF) for the double exponential distribution reads where the quantities and 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 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 grid of values and, for each pair, we used the NR scheme for classifying a set of initial conditions, on the configuration plane and in particular inside the squared region .
In Figure 6(a), we present the evolution of the average number of iterations , needed by the NR method for providing the coordinates of the equilibria with the desired accuracy. Figures 6(b) and 6(c) depict the distributions of the location parameter and the diversity 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 the fact that the numerical values of and are very close: ). Additionally, from the distribution of the diversity , shown in Figure 6(d), we can conclude that the probability histograms are very well organized around the average value , since in most of the cases the numerical value of the diversity is relatively low . Finally, Figure 6(d), we show how the differential entropy, defined as , evolves as a function of the values . It is seen that both quantities and have very similar parametric evolutions. If we take into consideration the combined information from all four panels of Figure 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.
(a)
(b)
(c)
(d)
Previously, in Figure 3, we have seen that there are certain regions on the -plane, 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, we shall compute the basin entropy [25, 27]. This modern tool indicates the fractal degree of a basin diagram by examining its topological properties. In Figure 7(a), we show the distribution of the numerical values of , 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.
(a)
(b)
Unfortunately, the transition between smooth and fractal boundaries cannot be determined by the basin entropy . The main reason for this drawback is that the basin entropy addresses the uncertainty to link a set of initial conditions to their 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 , 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 that if the basin boundary is fractal or not, by using the so-called “log 2 criterion,” with the sufficient condition, if , then the boundary is certainly fractal. The distribution of the values of , as a function of , is given in Figure 7(b). We see that when eight points of equilibrium exist, the basin boundaries on the -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 and the fractal dimension 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 Figure 8, we present the dependence of the uncertainty dimension 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 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 with respect to that of the basin entropy . This similarity can be explained by considering that these two computer-based analysis techniques are grounded on box-counting methodologies.
5. 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 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 -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, we needed, per grid, roughly 3 minutes using a Quad-Core i7 4.0 GHz CPU. All the plots of the paper have been developed using the software Mathematica [35].
Data Availability
The data underlying this article will be shared by the corresponding author upon reasonable request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was partially supported by COLCIENCIAS (Colombia) (Grant 8863) and by Universidad de los Llanos.