Solution (source code)

= Solution

The first <partial derivatives> are
$$
f_x=\frac1y-\frac{y}{x^2}-\frac{2(x-y)}{a^2},\qquad
f_y=\frac1x-\frac{x}{y^2}+\frac{2(x-y)}{a^2}.
$$
Both vanish at $(\lambda,\lambda)$, so every positive diagonal point is a <stationary point>. The <Hessian matrix> there is
$$
\boxed{H=2\left(\frac1{\lambda^2}-\frac1{a^2}\right)
\begin{pmatrix}1&-1\\-1&1\end{pmatrix}.}
$$
Its <eigenvectors> $(1,1)$ and $(1,-1)$ have respective <eigenvalues>
$$
\boxed{0,\qquad4\left(\frac1{\lambda^2}-\frac1{a^2}\right).}
$$
The zero <eigenvalue> reflects the entire line of <nonisolated stationary points>: $f(\lambda,\lambda)=2$. In fact $f-2=(x-y)^2[1/(xy)-1/a^2]$, showing a non-strict <local minimum> when $\lambda<|a|$, a non-strict <local maximum> when $\lambda>|a|$, and values on both sides of two near $\lambda=|a|$. At that last point the whole <Hessian matrix> vanishes, so its signs alone cannot classify the point.