I am able to accurately calculate the eigenvalues in perturbation theory, as follows:

I perform a similarity transformation on the matrix $M$, given by $M'=S^{-1}MS$ with
<sub>
$$S=\frac{1}{2m}\left(
\begin{array}{cccc}
 0 & 2 m & -\sqrt{m^2+4 n^2}+m-2 n & \sqrt{m^2+4 n^2}+m-2 n \\
 0 & 0 & -\sqrt{m^2+4 n^2}+m-2 n & \sqrt{m^2+4 n^2}+m-2 n \\
 2 m & 0 & 2 m & 2 m \\
 0 & 0 & 2 m & 2 m \\
\end{array}
\right).$$
</sub>
The matrix $M'$ has the same eigenvalues as $M$ and is given by
$M'=M_0+M_1$, with
<sub>
$$M_0=\left(
\begin{array}{cccc}
 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 \\
 0 & 0 & \frac{1}{2} \left(-\sqrt{m^2+4 n^2}+m+2 n\right) & 0 \\
 0 & 0 & 0 & \frac{1}{2} \left(\sqrt{m^2+4 n^2}+m+2 n\right) \\
\end{array}
\right),$$
$$M_1=\frac{1}{\sqrt{m^2+4 n^2}}\left(
\begin{array}{cccc}
 -2 \sqrt{m^2+4 n^2} & \sqrt{m^2+4 n^2} & -\frac{\sqrt{m^2+4 n^2} \left(\sqrt{m^2+4 n^2}-m+2 n\right)}{2 m} & \frac{\sqrt{m^2+4 n^2} \left(\sqrt{m^2+4 n^2}+m-2 n\right)}{2 m} \\
 \sqrt{m^2+4 n^2} & -\sqrt{m^2+4 n^2} & \sqrt{m^2+4 n^2} & \sqrt{m^2+4 n^2} \\
 \sqrt{m^2+4 n^2}-2 n & \frac{1}{2} \left(\sqrt{m^2+4 n^2}-m-2 n\right) & -\frac{3}{2} \sqrt{m^2+4 n^2}-\frac{m}{2}+n & -\frac{1}{2} \sqrt{m^2+4 n^2}-\frac{m}{2}+n \\
 \sqrt{m^2+4 n^2}+2 n & \frac{1}{2} \left(\sqrt{m^2+4 n^2}+m+2 n\right) & \frac{1}{2} \left(-\sqrt{m^2+4 n^2}+m-2 n\right) & \frac{1}{2} \left(-3 \sqrt{m^2+4 n^2}+m-2 n\right) \\
\end{array}
\right).$$
</sub>

Now we note that for large $m,n$ of order $1/\epsilon$, the matrix $M_0$ is of order $1/\epsilon$ while the matrix $M_1$ is of order unity. We can therefore perform perturbation theory in the small parameter $\epsilon$. 

Let me first look at the positive eigenvalues. To zeroth order these are given by $\beta_\pm=\frac{1}{2} \left(\pm \sqrt{m^2+4 n^2}+m+2 n\right)$. To first order these are the $(3,3)$ and $(4,4)$ diagonal elements of $M_0+M_1$, which gives the two positive eigenvalues
$$\beta_\pm=\frac{1}{2} \left(m+2 n\pm\sqrt{m^2+4 n^2}\right)\pm\frac{m-2 n
\mp 3 \sqrt{m^2+4 n^2}}{2 \sqrt{m^2+4 n^2}}.$$
This expression for the positive eigenvalues is quite accurate, see the plot for $m=n$, where the exact expressions (red and green lines) are almost indistinguishable from the perturbative result (blue and orange).

<IMG SRC="https://ilorentz.org/beenakker/MO/betaeigenavalues.png" WIDTH="400"/>

Now for the negative eigenvalues. To first order these are given by the eigenvalues of the $2\times 2$ upper-left block of $M_0+M_1$, so $\alpha_1,\alpha_2$ are the eigenvalues of the matrix ${{-2\;1}\choose{1\;-1}}$, given by
$$\alpha_1=\frac{1}{2} \left(-\sqrt{5}-3\right)=-2.618,\;\;\alpha_2=\frac{1}{2} \left(\sqrt{5}-3\right)=-0.382.$$

Here is the plot for $m=n$, again red and green is exact, blue and orange the perturbative result.

<IMG SRC="https://ilorentz.org/beenakker/MO/betaeigenvalues.png" WIDTH="400"/>


So to first order in $\epsilon$ one has $\alpha_1+\alpha_2=-3$. The challenge is to prove that higher order corrections are negative. Incidentally, the restriction to integer $m,n$ does not seem to play any role (in the plot $m=n$ is varied continuously).