Some notes on the Henon-Heiles Hamiltonian system

Posted on January 23, 2021 by mAnasa-taraMgiNI

Anyone familiar with dynamical systems knows of the Henon-Heiles (HH) system. What we are presenting here is well-known stuff about which reams of material have been written. However, we offer certain tricks for visualizing this system that make it easy for lay readers with just a high school knowledge of mathematics to play with. The HH system was discovered by the French astronomer Henon and his colleague Heiles when they were studying the motion of stars in the galaxy under the influence of the gravity of the total matter in the galaxy. The true astronomical significance of these equations outside the scope of the current discussion. Our own original interest in this problem was primarily from the perspective experimental mathematics (“play physics”), starting as an extension to our interest in defining ovals by means of ordinary differential equations (ODEs). The system defined by Henon and Heiles considers the motion of a body in 2 dimensional Euclidean space, i.e. a fixed plane. The phase space describing the motion is thus defined by the variables (x, y, p_x, p_y), where x, y describe the position in two dimensions and p_x, p_y describes the momentum in the two directions. Given that the momenta are p_x= mx'; p_y=my', for a body of unit mass the momenta become the derivative of the position variables with respect to time (x'=\tfrac{dx}{dt}, y'=\tfrac{dy}{dt}). Henon and Heiles considered a potential described by the equation:

V(x,y) = \dfrac{x^2+y^2}{2}+x^2y-\dfrac{y^3}{3} \kern 3em \cdots \S 1

The potential energy of a simple harmonic oscillator in the x direction is V(x)=\tfrac{kx^2}{2}. By taking a unit force constant k we see that the terms \tfrac{x^2+y^2}{2} in \S 1 represent two orthogonal simple harmonic oscillators. The further nonlinear term, x^2y-\tfrac{y^3}{3}, in \S 1 is a perturbation that couples these oscillators. This potential takes the form of a cubic hyperboloid-paraboloid and is visualized in Figure 1.

Figure 1.

The kinetic energy of the body is given by T=\tfrac{mv^2}{2}; where v is the velocity of the body. Thus, for the above-defined HH system we get T=\tfrac{x'^2+y'^2}{2}. The Hamiltonian of a system, which represents its total energy, is given by H=T+V. Since this is an energy conserving system, its total energy is equal to a scalar constant E, i.e. the energy level of the system. Thus, for the HH system we get:

H=\dfrac{x'^2+y'^2+x^2+y^2}{2}+x^2y-\dfrac{y^3}{3}= E \kern 3em \cdots \S 2

If we section the 3D curve V(x,y) by planes corresponding to different energy levels z=E, we get the equipotential curves within which the x,y would lie for a given energy level (Figure 2). We observe that if E=\tfrac{1}{6}, the equipotential curve becomes 3 intersecting lines that form an equilateral triangle defined by the equilibrium points (-\tfrac{\sqrt{3}}{2}, -\tfrac{1}{2}); (0,1); (\tfrac{\sqrt{3}}{2}, -\tfrac{1}{2})). Within this equilateral triangle, the body exhibits bounded motion. Thus, for all E<\tfrac{1}{6} we get bounded trajectories in the x-y plane. As E becomes smaller the central equipotential boundary tends towards a circle and degenerates to a point at 0. However, for E>\tfrac{1}{6} we get curves that are open; hence, at these energy levels the body can escape to infinity via the open lanes. Thus, there is a clearly defined escape energy level for this system, E=\tfrac{1}{6}.

Figure 2. The energy levels correspond to E=\tfrac{1}{32}, \tfrac{1}{16}, \tfrac{1}{12},\tfrac{1}{8},\tfrac{1}{6}, \tfrac{1}{4}, \tfrac{1}{3}, 1, 2, 3, 4

To study the trajectories under this system we first obtain the equations for the force acting on a body of unit mass (acceleration) in each direction from the above potential by taking the negative partial derivative with respect to each positional variable:

x''= -\dfrac{\partial V(x,y)}{\partial x} = -x-2xy

y''= -\dfrac{\partial V(x,y)}{\partial y} = -y -x^2+y^2

From the above can now get a system of ODEs thus:

x'=p_x
y'=p_y
p_x'=x''= -x-2xy
p_y'=y''=-y -x^2+y^2 \kern 3em \cdots \S 3

The solutions to this system \S 3 yield a curve in the 4-dimensional phase-space (x,y, p_x, p_y). To solve \S 3 we first need to obtain some initial conditions for a given energy level E using the Hamiltonian \S 2. We do that by setting x_0=0. We then choose some values of y_0, p_{y0}=y_0'. From those we can calculate p_{x0}=x_0' thus:

p_{x0}= \sqrt{2E-p_y^2-y^2+\dfrac{2y^3}{3}}

One can see that this places a constraint on the allowed y_0, p_{y0} — they have to be chosen such that p_{x0} is real. Once we have these initial conditions we can solve the above ODEs with efficient LSODA solver written by Alan Hindmarsh and Linda Petzold or you can write your own solver by the method of Runge and Kutta as we did in our youth. Initial results below are shown using the LSODA solver. However, we will see below that we can also obtain solutions without using traditional ODE solutions. Figure 4 shows an example of solution for the energy level E=\tfrac{1}{8} and initial conditions x_0=0; y_0= 0.1, y_0'= 0.14 in the 3D space defined by x, y, y'

Figure 3.

To get a better understanding of its behavior, we can visualize the solution in several other ways Figure 4. First, we can simply look at the way x, y change with time (first 2 top left panels of Figure 4). As expected, x(t), y(t) would be oscillatory functions that cannot be defined using any elementary functions. We can also examine the positional trajectory of the body in its plane of motion by plotting x, y (top right panel of Figure 4). From the equipotential curves defined above from \S 1, we can see that this trajectory would be bounded by the closed loop of the curve defined by the equation (shown in blue):

\dfrac{x^2+y^2}{2}+x^2y-\dfrac{y^3}{3}=E

We can also plot y, y' (bottom left panel of Figure 4) which shows how momentum changes with the position in the y direction. This curve will be bounded by a special oval (shown in blue) that is determined by letting x=0; x'=0 in the Hamiltonian \S 2. This gives us a cubic curve defined by the equation (in standard x-y coordinates, not the (x, y) of the phase space of the solutions of \S 3):

\dfrac{x^2+y^2}{2}-\dfrac{x^3}{3}=E \kern 3em \cdots \S 4

The closed loop of the cubic \S 4 is the bounding oval, which was what got us first interested in the HH system in the 16th year of our life.

Figure 4.

Finally, the bottom right panel of Figure 4 shows the Poincare section that records the points where the curve shown in Figure 3 pierces the plane x=0 (See below for further discussions). It is obvious that these are a subset of the y, y' plot and will thus be bounded by the same oval \S 4.

The way to compute the Poincare section is to search the x values of solution for cases where the sign of x_n and x_{n+1} are different. Such successive points will bound the segments of the curve that pierce the plane x=0. Given that our steps for numerical integration are small, we can calculate the corresponding values of y, y' using linear interpolation: y=\tfrac{y_n+y_{n+1}}{2}; \; y'=\tfrac{y_n+y_{n+1}}{2}. Plotting the thus calculated y, y' will give us the Poincare sections for a given starting point. Now, we can also calculate the solutions for above system \S 3 without solving the ODEs by converting it into a discrete difference equation. These difference equations have a step parameter \epsilon, which if kept small can yield solutions equivalent to that obtained by solving the ODEs. The system of difference equations goes thus:

p_{xn+1}= p_{xn}+\epsilon (-x_n-2x_ny_n)
p_{yn+1}= p_{yn} + \epsilon (-y_n+y_n^2-x_n^2)
x_{n+1}= x_n+\epsilon p_{xn+1}
y_{n+1}= y_n+\epsilon p_{yn+1}

Figure 5.

We empirically determined that by setting \epsilon =0.02 we can get results similar to the solution of the ODEs with time steps of 0.01. This provides us an easy mechanism, with somewhat higher speed than the ODE solver, to obtain equivalent solutions for the HH system. This in turn allows us to explore the Poincare sections for different initial values at a much higher density. Figure 5 shows one such exploration of Poincare sections for the energy level E=0.128 with 100 different initial conditions, each plotted in a different color. The result is a beautiful oval with an inner decoration by a strange attractor reminiscent of one of the ovoids produced for the Russian royalty. The attractor shows clear preferred regions for the intersections of certain orbits and regions where the intersections are chaotically distributed. To better understand the relationship between the structure of the Poincare sections and the form of the orbits on the x-y plane we take the case of E=\tfrac{1}{8} and examine 12 initial points chosen from different regions of the Poincare sections, i.e. defined y_0, y_0', with x_0=0 (Figure 6).

Figure 6.

The trajectories of these initial points on the x-y plane are plotted in Figure 7. Towards the narrow end of the bounding oval we have an oval exclusion zone and the towards the broad end of the oval we have a candra-bindu (crescent and dot) clearing zone. The initial point 1 lies at the center of the narrow end oval clearing. This initial point and the center of the crescent clearing at the broad end (not shown) yield a trajectory with a single loop with 3 apexes (top row, leftmost plot of Figure 7). The next trajectory (top row, next plot moving left to right) is a straight line at a 45^\circ incline and corresponds to the center of the two “eyes” of the Poincare section (point 2 in figure 6 is one of these eyes).

Figure 7. The trajectories of the points corresponding to Figure 6 in left to right in 3 rows from top to bottom.

The basis of these trajectories can be understood from the plots of the functions x(t); y(t) (Figure 8; for every point in Figure 6 and trajectory in Figure 7 the corresponding x(t); y(t) are shown one below the other from top to bottom in 2 columns). The first two trajectories result from oscillations where both x(t) and y(t) have period 1 — they show the same repeating pattern after one oscillation. Thus, these two cases can be said to be in a 1:1 resonance. In the second case, they are additionally in phase, i.e. the crest and trough at the same time.

Point 3 samples the center one of the four “islands” which surround the above-mentioned “eyes” of the Poincare section. Each of the island-centers results in a trajectory like the 3rd plot (Figure 7, top row). Point 4 samples one of the small crescents in the vicinity of the oval exclusion zone around point 1 and results in the trajectory seen in plot 4 of Figure 7. These two trajectories result from x(t); y(t) where both have a periodicity of 4, i.e. a 4:4 resonance. Of the two the trajectory 3 arises from a case where in addition to 4:4 resonance the two oscillators are also in phase.

Point 5 (Figure 6) and its corresponding trajectory (Figure 7) corresponds to two period 5 oscillators in a 5:5 resonance (Figure 8). Such 5:5 resonance oscillators are a pervasive feature of the HH system at this energy level and correspond to the 5 islands of exclusion around the oval exclusion around point 1, the center of the bindu and the two exclusion zones flanking either tip of the crescent.

Figure 8

Point 6 corresponds to a trajectory arising from a 8:9 resonance; point 7 evolves into a more complex 5:25 resonance; the trajectory of point 8 simulates a 3D ribbon and arises from the even more complex 11:37 resonance.

The trajectories arising from points 9, 10 and 11 exhibit what might be termed quasiperiodic behavior. In the case of the evolution of point 9, x(t) has a quasiperiod of 4, i.e., it has a similar pattern repetition after every 4 oscillation but the successive repeats are not identical but change slightly over time. y(t), on the contrary, has a strict period of 1. In the evolution of point 10, x(t) has a quasiperiod of 5 which is overlayed on a nearly regular higher period of 15. These two points are representative of the evolution of the points in the zones close the bounding oval on its narrow side. One may note that the evolution of point 11 is like a “broadband” version of the 5:5 resonance points. Keeping with this, x(t) has a strict period of 5, whereas y(t) has a quasiperiod of 5 with higher-order repeat patterns of multiples of 5.

Finally, the evolution of point 12 is chaotic, i.e. the oscillations have no discernible period. The irregularity is marked in y(t) but is x(t) it manifests more slowly over time. These chaotic trajectories form the bulk of the central uniform distribution of points in the Poincare section. The appearance of chaos can be seen as the limit of the trajectories with increasingly complex or longer period resonances. The quasiperiodic orbits with a nearly regular short period internal repeat structure might be seen as lying at the edge of long periods and true periodicity. In terms of energy levels, chaos starts appearing in the central regions close to E=\tfrac{1}{10} and by E=\tfrac{1}{8} constitutes the bulk of the internal structure of the Poincare section with internal islands of periodicity and quasiperiodicity in the anterior periphery of the oval. By the limiting E=\tfrac{1}{6} nearly all of the trajectories become chaotic.

In conclusion, the HH system qualitatively shows all the typical forms of oscillatory behaviors observed in natural systems (e.g. variable star light curves, weather patterns, population dynamics and far-from equilibrium oscillatory chemical reactions): periodicity with different resonances, quasiperiodicity and chaos. It thus provides a good example how any system whose phase space is defined by even simple ODEs with non-linear terms can exhibit the behavioral diversity characteristic of natural systems.

This entry was posted in Scientific ramblings and tagged , , , , , , , , , , , , , . Bookmark the permalink.