This is not yet an answer but some progress: I am able to accurately calculate the two positive eigenvalues, 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$. To first order the eigenvalues of $M$ are the diagonal elements of $M_0+M_1$, which gives the two negative eigenvalues
$\alpha_1=-2$, $\alpha_2=-1$, and the two positive eigenvalues $\beta_+$ and $\beta_-$ given by 
$$\beta_\pm=\frac{1}{2} \left(\sqrt{m^2+4 n^2}+m+2 n\right)\pm\frac{m-2 n-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"/>

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).