Henon Heiles

Published

October 30, 2026

Also available as a slide deck.

The Henon–Heiles describes a celestial system of a star inside of a effective potential of a galaxy. It is a well understood and results can be compared to existing results.

Hénon-Heiles potential

\[\Epot(x,y)=\frac{m\omega^2}{2}(x^2+y^2)+\alpha(x^2y-\frac{y^3}{3})\]

Hamiltonian & equations of motion

\[\begin{align} H(x,y,p_x,p_y) &= \Ekin(p_x, p_y) + \Epot(x,y) \\ =&\frac{(p_x^2+p_y^2)}{2m} + \frac{m\omega^2}{2}(x^2+y^2)+\alpha(x^2y-\frac{y^3}{3}) \\ \overset{\alpha=1, m=1, \omega=1}{\longrightarrow}\, &\frac{(p_x^2+p_y^2)}{2} + \frac{(x^2+y^2)}{2}+(x^2y-\frac{y^3}{3}). \end{align} \]

We obtain the equations of motion \((\dot{x_i}\equiv\tfrac{\partial x_i}{\partial t})\) by

\[ \dot{q_i} =\frac{\partial H}{\partial p_i} = \frac{\partial\Ekin}{\partial p_i} \]

\[ \dot{p_i} =-\frac{\partial H}{\partial q_i} = -\frac{\partial \Epot}{\partial q_i} \\ \]

canonical conjugate variables \((p_i,q_i)_{i\in\mathbb{N}}\) labeled by \(i\in (1,2 \dots, n)\).
In our case simply \(q_i\in\{x,y\}\) and \(p_i\in\{p_x, p_y\}\) span the phase space
phase space: \(\textbf{x}=(\textbf{p}, \textbf{q}) \in \mathbb{R}^N\), \(N=4\).

The equations of motion are then derived by \(\dot{q_i} =\tfrac{\partial H}{\partial p_i}\) and \(\dot{p_i} = -\frac{\partial \Epot}{\partial q_i}\) resulting in

\[\begin{align} x&=\frac{p_x}{m} &\longrightarrow & &p_x \\ y&= \frac{{p_y}}{m} &\longrightarrow & &p_y\\ p_x&= -m\omega^2\,x - 2\alpha\, xy &\longrightarrow & & -x-2xy \\ p_y&= -m\omega^2 y - \alpha (x^2-y^2) &\longrightarrow & &-y- x^2+y^2 \end{align}\]

It may be helpful to define the following:

\(H(\textbf{x},t)= H(\textbf{p},\textbf{q},t)\equiv H(p_x,p_y,x,y)\),
\(\Ekin(\textbf{p})\equiv\Ekin(p_x, p_y)=\Ekin(p_x)+\Ekin(p_y)\),
\(\Epot(\textbf{q})\equiv\Epot(x,y)\) and finally
\(\Epot(y)=\Epot(0,y)\) or \(\Epot(x)=\Epot(x,0)\)

Pioncré section

To understand our 4-dimensional phase space:

the system is bound by Energy \(E\) available \[H(\textbf{x})=E=\Ekin(\textbf{p})+\Epot(\textbf{q})\]
Due to \(\Ekin(p_i)\geq 0\) \[\Rightarrow \;E\geq\Ekin(p_y)+\Epot(\textbf{q})\] with equality if and only if \(p_x=0\)
to check for periodic orbits we create surface of section \(\sufsec\in\mathbb{R}^2\) at \(x=0, p_x\geq0\), \[E\geq\Ekin(p_y)+\Epot(y)\]
find the boundary of the possible phase space points passing \(\sufsec\), i.e. \[\{(y,p_y) | E=\Ekin(p_y)+\Epot(y)\}\]

Implementation

The resulting equation of motion in julia code:

function equations!(du, u, p, t)
    a, m, w = p 
    x, y, px, py = u
    du[1] =  px
    du[2] =  py
    du[3] = -m*w^2*x - 2a*x*y
    du[4] = -m*w^2*y - a*(x^2 - y^2)
end

Here the potential and kinetic energy are defined:

potential(x, y, p) = p[2]*p[3]^2/2 *(x^2+y^2) + p[1]*(x^2*y - y^3/3)

kinetic(px, py, p) = (px^2 + py^2)/(2*p[2])

hamiltonian(x, y, px, py, p) = kinetic(px, py, p) + potential(x,y, p)

The equations of motion are solved by use of the package(Datseris and Parlitz 2022). The data for the surface of section is gathered in runtime of the solver. This is achieved by creating a ‘callbackfuntion’, conditioned on \(x=0, p_x\geq0\), appending \((y,p_y)\) to their respective lists:

function section_callback(section_q, section_p)
    condition(u, t, integrator) = u[1]
    affect!(integrator) = begin
        if integrator.u[3] >= 0
            push!(section_q, integrator.u[2])
            push!(section_p, integrator.u[4])
        end
    end
    ContinuousCallback(condition, affect!)
end

The trajectory gets solved as followed

function solve_trajectory(u0, p, tspan, dt; callback=nothing)
    prob = ODEProblem(equations!, u0, tspan, p)
    solve(prob, Vern9(), abstol=1e-10, reltol=1e-10, saveat=dt, callback=callback)
end

Results

Due to \(E\geq \Ekin(p_y)+\Epot(y)\) we can obtain a boundary of our plot for our surface of section, which correspond to \(x=0, p_x=0\).

solve \(E-(\Ekin(p_y)+\Epot(y))=0\Rightarrow (p_y,y)\)?
\(p_y=0\) to obtain 3rd-order polynomial: \(E-\Epot(y)=0\Rightarrow (y_\text{min}, y_\text{max})\)
numerically determine \(E-(\Ekin(p_y)+\Epot(y))=0\Rightarrow (p_y,y)\) for \(y\in[y_\text{min}, y_\text{max}]\)

The results presented consistently use the parameters

\(m=1\), \(\omega=1\), \(\alpha=1\) and saves at timesteps \(dt=0.5\)
the goal was to identify various periodic trajectories manually for \(E=0.1127\)

References

Datseris, George, and Ulrich Parlitz. 2022. Nonlinear Dynamics: A Concise Introduction Interlaced with Code. Springer Nature. https://doi.org/10.1007/978-3-030-91032-7.

Hénon-Heiles potential

Hamiltonian & equations of motion

Pioncré section

Implementation

Results

the goal was to identify various periodic trajectories manually for \(E=0.1127\)

References