# Discovery of Maximally Super Integrable Quantum Mechanical Systems

## Mathematical Context
Suppose one is given a Hamiltonian operator
$$
H=-\partial_x^2+V(x)
$$
on $\mathbb{R}^n$ with real-valued potential $V(x): \mathbb{R}^n \rightarrow \mathbb{R}$. An "integral" of motion is given by an operator $K_i$ on $\mathbb{R}^n$ such that the matrix commutator
$$
[H,K_i]=0
$$
holds for all functions $\psi\in L^2(\mathbb{R}^n)$. These integrals may be algebraically dependent, so to ensure that we have genuinely unique integrals, we demand that each integral does not commute with one another, i.e.,
$$
[K_i,K_j]\neq0\qquad \forall \ i\neq j.
$$
If there are $n$ integrals (including the Hamiltonian $H$ which is itself an integral of the motion), we say that the system is integrable. If we have $2n-1$ we say that the system is maximally super-integrable.

As a concrete example of a maximally super-integrable system, consider the two dimensional family of Calogero models, parametrized by $c,\omega\in\mathbb{R},$ and given by the Hamiltonian

$$
H=-\partial_{q_1}^2-\partial_{q_2}^2+\omega^2\left(q_1^2+q_2^2\right)+\frac{2 c^2}{\left(q_1-q_2\right)^2} \\
$$
We observe that the following integrals $K_1$ and $K_2$
$$
\begin{aligned}
& K_1=\left(q_2 \partial_{q_1}-q_1 \partial_{q_2}\right)^2-\frac{4 c^2 q_1 q_2}{\left(q_1-q_2\right)^2} \\
& K_2=\partial_{q_1} \partial_{q_2}-\omega^2 q_1 q_2+\frac{c^2}{\left(q_1-q_2\right)^2}
\end{aligned}
$$
satisfy the above conditions for algebraically independent integrals of motion.

The point of this research project is to algorithmically discover integrals for systems that are so high-dimensional that even computer algebra systems like Maple and Mathematica are rendered useless.

## Problem Formulation and Numerical Framework

We relax the exact calculation of commutation relationships among operators to computations that are handled entirely via numerical optimization and random sampling techniques. To help keep our framework concrete, we will specialize to the two-dimensional Calogero model. Generalizing the method we will describe below is a straightforward task.

### Operator Hypothesis
Let us hypothesize that the first and second integrals are of the form
$$
\tilde{Q}=\sum_i\sum_j\sum_k\sum_l A_{ijkl} q_1^iq_2^j\partial_{q_1}^k\partial_{q_2}^l+\frac{\sum_i\sum_jB_{ij}q_1^iq_2^j}{\sum_i\sum_jC_{ij}q_1^iq_2^j}
$$
The minimal number of parameters introduced so that this hypothesis is wide enough to enclose the actual integrals $K_1$ and $K_2$ is $2\times(3^4+3^2+3^2)=198.$ This problem is already fairly taxing. In any case, we proceed forward and use operator hypotheses of the form $\tilde{Q}$ to discover the hypothesized integrals $\tilde{K}_1$ and $\tilde{K}_2$ numerically.

### Sampling over the Function Space
To construct the loss function that models the discovery of integrals, we must also mimic the procedure of evaluating the commuation relations pointwise over function spaces. to this end, knowledge of the diagonalization of the Hamiltonian is of use. The eigenfunctions of the two-dimensional Calogero model are given by

$$
\varphi_m\left(q_1, q_2\right)=\left(q_1-q_2\right)^{\nu} e^{\frac{\omega}{2}\left(q_1^2+q_2^2\right.)} P_m\left(q_1, q_2\right)
$$

where $\nu=c^2(c^2-1)\geq-1/4$ ensures that each $\varphi_m\in L^2(\mathbb{R}^2)$ and $P_m$ is an $m^{\rm th}$ degree orthogonal polynomial whose explicit form is easy to determine via induction. To evaluate our commutators, we construct random functions that live in the span of $\{\varphi_m\}_{m=0}^{\infty}$. That is, use
$$
\psi_j=\sum_ma_m\varphi_m(q_1,q_2)
$$
where each $a_m$ is sampled from the uniform distribution $U[-1,1]$. Of course, we must truncate the sum. Twenty basis functions here should suffice for a thorough representation of the function space.

### The Loss Function
The first contribution to the loss function is the one that models the discovery of an integral which commutes with the Hamiltonian:
$$
J_{\rm C}=\sum_i^{2n-2}\sum_j\int_{\mathbb{R^2}}dq\left|[H,\tilde{K}_i]\psi_j\right|^2
$$
To model algebraic independence, we introduce the penalization
$$
J_{\rm A}=\sum_{i\neq k}\sum_j\frac{1}{\int_{\mathbb{R^2}}dq\left|[\tilde{K}_i,\tilde{K}_k]\psi_j\right|^2}
$$
If this loss is large, it means that the integrals are close to being algebraically dependent.

To promote a parsimonious discovery of superintegrability, we introduce a sparsification on the parameters via a thresholded regularization $J_{\rm S}$, please see this paper for more details:
https://arxiv.org/abs/2503.00645
This loss is known to promote concentration of the "mass" of the weight vector to a few important locations. This aids in interpretability of the numerical results.

### The Proposed Optimization Problem
The optimization problem we thus seek to solve is
$$
\min_{\xi\in\mathbb{R}^{198}}(1-r)J_{\rm C}[\xi]J_{\rm A}[\xi]+rJ_{\rm S}[\xi],
$$
where $r$ is our sparsifying regularization parameter and $\xi$ is the vectorization of the joint tensors $A,\ B,$ and $C$. The choice of this parameter value is ad hoc and dependent on computer experimentation. To solve this optimization problem, we propose to use stochastic gradient descent, with ADAM, and gradients computed via automatic differentiation. Much of this optimization machinery is incredibly standard and we will leverage the vast number of tools that PyTorch offers to this end.

## First Computational Results

I had quite a humbling experiencing with my first computations on this project. I somehow managed to overlook that **the Hamiltonian is singular.** This is turning out to be disastrous.

I set $\omega=2$ and $c=1$. As test data, I tried to evaluate the Hamiltonian acting on
$$
f(q_1,q_2)=(q_1-q_2)^{\nu}e^{-\omega/2(q_1^2+q_2^2)}(\sin(6q_1)+\cos(3q_2))
$$
You see, the Hamiltonian is singular along the line $q_1=q_2$

I need to figure out how I could possibly make sense of this when I evaluate the objective function for this problem. Anyway, here's the MATLAB code I used

%spatial grid
l=1;N=2^8;
x=l*(2*pi/N)*(-N/2:N/2-1)';
dx=x(2)-x(1);

%Build the wavenumber
dk=2*pi/(N*dx);
k=[0:dk:N/2*dk,-(N/2-1)*dk:dk:-dk]';

% q grid
[q1,q2]=meshgrid(x);

% system parameters
w=2;c=1;

% test data
F=getEigenFun(q1,q2,c,w);

% This takes partial derivatives
D=@(u,which,m)TakePartialDerivative(u,k,which,m);

% the Hamiltonian
H=@(u)-D(u,1,2)-D(u,2,2)+w^2*(q1.^2+q2.^2)+2*c^2./(q1-q2).^2;

figure
pcolor(q1,q2,real(F)), shading interp, colorbar
xlabel('$q_1$'),ylabel('$q_2$'), title('$f(q_1,q_2)$')
figure
pcolor(q1,q2,real(H(F))), shading interp, colorbar
xlabel('$q_1$'),ylabel('$q_2$'), title('$Hf(q_1,q_2)$'), clim([-30 80])

function F=getEigenFun(q1,q2,c,w)
nu=c^2*(c^2-1);
peel=(q1-q2).^nu.*exp(-w/2*(q1.^2+q2.^2));
F=peel.*(sin(6*q1)+cos(3*q2));
end

function Fq=TakePartialDerivative(F,k,which,m)
% computes spectral derivatives
D=@(u,m)TakeSpectralDerivative(k,u,m);
Fq=0*F;
N=length(k);
if which==1
    for i=1:N
        % this is a row derivative
        Fq(i,:)=D(F(i,:).',m).';
    end

else
    for i=1:N
        % this is a column derivative
        Fq(:,i)=D(F(:,i),m);
    end
end
end

function up=TakeSpectralDerivative(k,u,m)
up=(ifft(((1i*k).^m).*fft(u)));
end

## A Strategy to Regularize

Of course, there's an obvious strategy I could use (at least it became clear to me after some coffee on March 26th). Let's say, anywhere I have a singularity that looks like

$$
\mathcal{O}_{\rm singular}=\frac{\mathcal{O}_{\rm just\ fine}(\vec{q},\nabla,\nabla \vec{q})}{|q_i-q_j|^p},
$$
where $\mathcal{O}_{\rm just\ fine}$ is an operator that will evaluate just fine, I could instead modify this so that instead, I evaluate the regularized operator
$$
\mathcal{O}_{\varepsilon}=\frac{\mathcal{O}_{\rm just\ fine}(\vec{q},\nabla,\nabla \vec{q})}{|q_i-q_j|^p+\varepsilon},
$$
I could then at least controllably evaluate the singularity of the operator. I found by experimenting, that

$\varepsilon=.01$ looks like this

$\varepsilon=.1$ looks like this

and $\varepsilon=.5$ looks like this

I can work with this! I should make sure that anywhere I have a singularity, I use the same parameter $\varepsilon$ to regularize. I will then have to experiment with this parameter to see how to balance this with the true evaluation of the singular operators when testing for superintegrability