# A Phase Space Study of the Nearly Integrable Henon-Heiles Model Using the Lax Spectrum

## Context
There is a general desire to understand, with great detail, the breakdown of integrability in Hamiltonian dynamical systems. People have been interested in this far before I was even born, so I feel fortunate that this interest has been passed to me.

<!-- Truly excellent work has recently been undertaken by Wei and Panos toward this end. Moreover, Yannis and his postdoc, who I only met just last night, have been spearheading studies to understand the limits to so-called "flow boxing." I only heard about flow boxing for the first time in March. Here's a no nonsense reference that I read to get me up to speed on what flow boxing is
https://arxiv.org/pdf/math/0305207 -->

Wei and Panos have used scientific machine learning to detect when integrability breaks down. Their approach is fairly structure agnostic, and hence, can be applied quite broadly to various dynamical systems of interest.

Being the pedestrian that I am, I chose to look into the Henon-Heiles (HH) system. The HH Hamiltonian is given by
$$
H=\frac{1}{2}\left(p_x^2+p_y^2\right)+\frac{1}{2}\left(A x^2+B y^2\right)+x^2 y+C y^3.
$$
Integrability is known to depend on the parameters $(A,B,C)$ defining the Hamiltonian. See this French publication for more details
https://www.math.univ-toulouse.fr/~gavrilov/publications/18.pdf

**Warning**: Do not completely trust this publication because it has some errors. Nevertheless, it's kind of the best reference I could find for the Henon-Heiles model. Anyway, we will focus on the case $A=B$ and $C=\frac{1}{3}$, in which case, a Lax pair is known:
\begin{equation}
L=\left(\begin{array}{cc}
L_1 & 0 \\
0 & L_2
\end{array}\right), \quad P=\left(\begin{array}{cc}
P_1 & 0 \\
0 & P_2
\end{array}\right)
\end{equation}
where
\begin{equation*}
\begin{aligned}
& L_1=\left(\begin{array}{cc}
\frac{p_y-p_x}{2} & y-x-k \\
\frac{2 k^2+k[3 A+2(y-x)]+3 A(y-x)+2(y-x)^2}{12} & \frac{p_y-p_x}{2}
\end{array}\right), \quad P_1=\left(\begin{array}{cc}
0 & 1 \\
-\frac{y-x}{3}-\frac{3 A+2 k}{12} & 0
\end{array}\right), \\
& L_2=\left(\begin{array}{cc}
\frac{p_y+p_x}{2} & y+x-k \\
\frac{2 k^2+k[3 A+2(y+x)]+3 A(y+x)+2(y+x)^2}{12} & -\frac{p_y+p_x}{2}
\end{array}\right),\quad P_2=\left(\begin{array}{cc}
0 & 1 \\
-\frac{y+x}{3}-\frac{3 A+2 k}{12} & 0
\end{array}\right),
&
\end{aligned}
\end{equation*}
and where $k$ is a free parameter. Notice that the momenta have a factor of two in this note that is missing from the paper by Gavrilov and others. That factor of two cost me an entire work day once, but I think it's probably time for me to stop crying about it now.

## Approach and Numerical Results
The Henon-Heiles system is truly a wonderful system to study. It is such a simple two degree of freedom system that offers a playground to benchmark our ideas, including the one I will demonstrate now. When $C=\frac{1}{3},$ the HH system has two conserved quantities: $CQ_1=H$ and
$$
CQ_2=\mathrm{Tr}\ L^4.
$$
For clarity, consider the Hamiltonian
$$
H=\frac{1}{2}\left(p_x^2+p_y^2\right)+\frac{1}{2}\left(x^2+y^2\right)+x^2 y+ \frac{1+\varepsilon}{3}y^3.
$$
As we vary $\varepsilon$ away from zero, we expect that the second conserved quantity should be destroyed. To model this destruction, let us consider the function
$$
\mathcal{C}_j(\varepsilon)=\int_{t_0}^{t_\rm f}\left(\frac{d}{dt}\mathrm{Tr}L^{2j}(x^{\varepsilon},y^{\varepsilon},p^{\varepsilon}_x,p^{\varepsilon}_y)\right)^2dt
$$
where the variables with superscripts are numerically integrated dynamics for a given value of $\varepsilon$, and where
$$
L=\left(\begin{array}{cc}
L_1 & 0 \\
0 & L_2
\end{array}\right)
$$
with
\begin{equation*}
L_1=\left(\begin{array}{cc}
\frac{p_y-p_x}{2} & y-x \\
\frac{y-x}{4}+\frac{(y-x)^2}{6} & \frac{p_x-p_y}{2}
\end{array}\right),\ \ L_2=\left(\begin{array}{cc}
\frac{p_y+p_x}{2} & y+x \\
\frac{y+x}{4}+\frac{(y+x)^2}{6} & -\frac{p_y+p_x}{2}
\end{array}\right).
\end{equation*}
The parameter $k$ was free, right? I set it equal to zero to keep things clean. Here's what $\log \mathcal{C}_1(\varepsilon)$ looks like from dynamics evolved from the initial condition $(x^{\varepsilon}(0),y^{\varepsilon}(0),p^{\varepsilon}_x(0),p^{\varepsilon}_y(0))=(.15,-.15,0,0)$ over the time domain $(t_0,t_{\rm f})=(0,1/2)$.

You would think I ripped this from the breakdown of KdV integrability that SILO found for us back in December (see Figure 5)
https://arxiv.org/abs/2503.00645
And yet, $\log\mathcal{C}(\varepsilon)$ is even smoother than that PDE example. It's important to note that this study was done for a single point in phase space, and doesn't inform us, as of yet, of anything related to the global structure of integrability

Now, this is what I observe if, instead, I compute $\mathcal{C}_2(\varepsilon)$ over $(x_0,y_0,0,0)$:
I observe a beautiful geometry emerging, perhaps indicating to us when an effective "one-body" description is no longer valid for a two-body problem. As expected, this seems to be valid and robustly persists in a wide neighborhood about the origin as nonintegrability is cranked up. Surprisingly, we see the emergence of algebraic curves where the second conserved quantity seems to remain valid. Please note that the scale is logarithmic, so -37 is about machine epsilon.

By the way, here is what $\mathcal{C}_1(\varepsilon)$ looks like integrated on $[0,1]$

#### Thoughts
* This picture looks different when the initial momenta are changed. This is a four dimensional phase space. I am not at the moment that motivated to look into this large dimensional space without more mathematical or physical understanding.
* I do not have a mathematical understanding of my computational results. This understanding is highly desirable to me (possibly us).
* I hope to have this effort here understood as complementary to the SciML approaches that are being undertaken by our study group. I wonder if you could use this HH case study as a guide towards refining your numerical results.
* The Lax pair is just so powerful, and this was fun to look into. I wonder what else we could say about other nearly integrable dynamical scenarios using the Lax pair as a proxy. Moreover, I believe this strengthens the notion of using SILO as a nonintegrability filter in our future endeavors.

## Code
Not the Cleanest Matlab Code, but it's here for reproducibility and transparency.

%% Test Study to Confirm Integrability
% perturbation values away from HH integrability
veps=-.1:.00001:.1;
J=0*veps;

for j=1:length(veps)

% Parameters
A       = 1;     % Coefficient for x^2 term
B       = 1;     % Coefficient for y^2 term
epsilon = 1/3+veps(j);   % Coefficient for the y^3 term

% Initial Conditions
% Choose initial conditions for x, y, p_x, p_y
x0  =  0.15;
y0  = -0.15;
px0 = 0;
py0 = 0;
z0  = [x0; y0; px0; py0];

% Time Span for Simulation
tspan = [0 2];

% get the Solution State
[t,z]=getState(A,B,epsilon,tspan,z0);

J(j)=getLoss(t,z);
end

plot(veps,log(J)), axis tight
xlabel('$\varepsilon$'),ylabel('$\log \mathcal{C}$')
grid on

%% FlowBox Study with GIF Output

% Domain Size
lxl = -1.02; lxr = 1.02;
lyl = -1.52; lyr = 0.3;
sx = 1/150; sy = sx;

% Initial Conditions
x0vec = lxl:sx:lxr;
y0vec = lyl:sy:lyr;
px0 = 0; py0 = 0;
N = length(x0vec);
M = length(y0vec);

% Time Span for Simulation
tspan = [0, 0.5];

% GIF Setup
gif_filename = 'HHflowbox.gif';
veps = 0.00000001;
numTrials = 50;

for k = 1:numTrials
    J = zeros(N, M) + 1e-16;
    epsilon = 1/3 + veps;

    for i = 1:N
        for j = 1:M
            % compute the loss at one phase space point
            z0 = [x0vec(i); y0vec(j); px0; py0];
            [t, z] = getState(A, B, epsilon, tspan, z0);
            J(i,j) = getLoss(t, z);
        end
    end

    % Plot
    pcolor(x0vec, y0vec, log(J)'), shading interp
    colorbar
    xlabel('$x_0$', 'Interpreter', 'latex')
    ylabel('$y_0$', 'Interpreter', 'latex')
    title(sprintf('$\\log \\mathcal{C}$ when $\\varepsilon = %.2e$', veps), ...
        'Interpreter', 'latex')
    colormap('jet')
    caxis([-40 15])
    drawnow

    % Capture frame and write to GIF
    frame = getframe(gcf);
    im = frame2im(frame);
    [imind, cm] = rgb2ind(im, 256);

    if k == 1
        imwrite(imind, cm, gif_filename, 'gif', 'Loopcount', inf, 'DelayTime', 0.3);
    else
        imwrite(imind, cm, gif_filename, 'gif', 'WriteMode', 'append', 'DelayTime', 0.3);
    end

    % Update epsilon
    veps = veps * 1.5;
end

%% Functions Used
function [t,z]=getState(A,B,epsilon,tspan,z0)
% Define the Hamiltonian ODE System
% The state vector is z = [x; y; p_x; p_y].
hamiltonianODE = @(t, z) [ z(3);                   % dx/dt = p_x
                           z(4);                   % dy/dt = p_y
                           -A*z(1) - 2*z(1)*z(2);   % dp_x/dt = -A*x - 2*x*y
                           -B*z(2) - z(1)^2 - 3*epsilon*z(2)^2]; % dp_y/dt

% Solve the ODE using ode45
[t, z] = ode45(hamiltonianODE, tspan, z0);
end

function J=getLoss(t,sol)

N=max(size(sol));
lambda=zeros(N,2);

% Compute the Traces of the Spectral operator over time
for k=1:N
    z=sol(k,:);
    L1=[z(4)/2-z(3)/2                           z(2)-z(1); 
        (3*(z(2)-z(1))+2*(z(2)-z(1))^2)/12  z(3)/2-z(4)/2];
    L2=[z(4)/2+z(3)/2                           z(2)+z(1); 
        (3*(z(2)+z(1))+2*(z(2)+z(1))^2)/12 -z(3)/2-z(4)/2];
    L=[L1 zeros(2);zeros(2) L2];

    % the second base conserved quantity we want to track
    lambda(k)=trace(L^4);
end

% this is the loss function
J=trapz(t,gradient(lambda,t).^2);
end