Asymmetric Hessians

I encountered a really counter-intuitive result this past week regarding the Hessians of a simple 2D Harmonic Oscillator:

$U(r(x_0, y_0, x_1, y_1), b) = (r-b)^2/2 =(\sqrt{(x_0-x_1)^2+(y_0-y_1)^2}-b)^2/2$

It is clear that the energies for $r+\delta$ and $r-\delta$ are the same, and the forces are identical in magnitude except opposite in direction.

However, while the internal Hessian $\dfrac{\partial^2 U}{\partial r^2} = 1$ and therefore positive semi-definite (PSD) every where w.r.t. $r$, this is not the case for the 4×4 Cartesian Hessian $\dfrac{\partial^2 U}{\partial x_i \partial x_j}$. In fact, the Cartesian Hessian is PSD only if $r \geq b$, that is, it has only positive eigenvalues or zero. We can actually analytically show that the eigenvalues of the Cartesian Hessian are $0, 0, 2$ and $\lambda$ where:

$\lambda= 2(1-\dfrac{b}{r^2})$

So if $r < b$ then the eigenvalue is negative (no longer PSD), and if $r>b$ then the eigenvalue is positive (PSD).

To see this for yourself in sympy, the code is:

import sympy as sp
from sympy.functions import Abs, sqrt
from IPython.display import display, Math
x0, x1, y0, y1, b = sp.symbols('x0 x1 y0 y1 b', real=True)
U = ((sqrt((x0-x1)**2+(y0-y1)**2)-b)**2)/2
sp.init_printing()
args = (x0, y0, x1, y1)
H = sp.Matrix(list(sp.derive_by_array(sp.derive_by_array(U, args), args)))
H = H.reshape(4,4)
H = sp.simplify(H)
H.eigenvals()