Eigenvalues of $\left[\begin{matrix}1&\epsilon&2\epsilon\\\epsilon&1&2\epsilon\\-2\epsilon&-2\epsilon&2\end{matrix} \right]$ as a series in $\epsilon$ Consider the matrix
$$
M := \left[ \begin{matrix} 1 & \epsilon & 2\epsilon \\ \epsilon & 1 & 2\epsilon \\ -2\epsilon & - 2\epsilon & 2 \end{matrix} \right]
$$
for some small parameter $0 < \epsilon \ll 1$. My question is how does one solve for the eigenvalues $\lambda$ of this matrix as a series in $\epsilon$, given that the zeroth-order (with $\epsilon = 0$) eigenvalues have degeneracy?
I am more interested in the procedure than anything else here. The eigenvalues are exactly given by
$$
\lambda \in \left\{ 1 - \epsilon, \frac{3 + \epsilon - \sqrt{ 1 - 2 \epsilon - 31 \epsilon^2}}{2} , \frac{3 + \epsilon + \sqrt{ 1 - 2 \epsilon - 31 \epsilon^2}}{2} \right\} \ ,
$$
which when expanded for $0  < \epsilon \ll 1$ gives respectively
$$
\lambda \in \left\{ 1 - \epsilon, 1 + \epsilon + 8 \epsilon^2 + \mathcal{O}(\epsilon^3)  , 2 -  8\epsilon^2 + \mathcal{O}(\epsilon^3) \right\} \ .
$$
How would one solve for these eigenvalues as a series if it was not possible to solve for the exact eigenvalues?
My attempt: Normally what I would do would be to assume I have a series
$$
\lambda = \lambda_0 + \epsilon \lambda_1 + \epsilon^2 \lambda_2 + \ldots
$$
and insert this series into the characteristic equation for the matrix which is here
$$
\lambda^3 - 4 \lambda^2 + (5 + 7 \epsilon^2) \lambda -2 - 6 \epsilon^2 + 8 \epsilon^3 = 0 \ .
$$
However this method yields nonsense because of the degeneracy of the first two eigenvalues (both being 1) in the limit that $\epsilon = 0$. To describe what happens: the characteristic equation tells me that $\lambda_0^3 - 4 \lambda_0^2 + 5 \lambda_0 - 2 = 0$ for $\epsilon =0$, which gives me that $\lambda_0 \in \{1,1,2\}$. However, at next order in $\epsilon$ in the characteristic equation tells me $5 \lambda_1 - 8 \lambda_0\lambda_1 + 3 \lambda_0^2 \lambda_1 = 0$, which when evaluated for the degenerate $\lambda_0 = 1$ just says that $0=0$ and so gives no information at all, and so I cannot proceed to find $\lambda_1$.
 A: This is really a question about the polynomial, not about the matrix. Consider a simple example $$(x-1)^2-a \epsilon^2 = 0$$ and expand: $x=1+x_1 \epsilon+ x_2 \epsilon^2 \cdots $. The first order in $\epsilon$ is indeed $0=0$. This does not mean you have "nonsense", it just means that $x_1$ is determined by the next order: $$x_1 ^2 \epsilon^2 -a \epsilon^2 =0$$
Similarly, for your more complicated example:
$$\lambda^3 - 4 \lambda^2 + (5 + 7 \epsilon^2) \lambda -2 - 6 \epsilon^2 + 8 \epsilon^3 = 0 \ $$
expanding to second order around $\lambda = 1 + \lambda_1 \epsilon +\cdots$ will give $1-(\lambda_1)^2 = 0$ as expected.
The reason is the double zero: there is a parabola $y=x^2$ that just touches the $x$ axis. An order $\epsilon$ perturbation in the constant $x^2 -\epsilon$ yields a $O(\sqrt \epsilon)$ in the root. Luckily for you, the perturbation is $x^2- \epsilon^2$, so you get an $O(\epsilon)$ change in the root position.
Edit
As an example of how things could have been even worse, consider the matrix $$\pmatrix{1&\epsilon\\a &1} $$ with the polynomial $$(1-\lambda)^2-a\epsilon$$ the lowest order term in the expansion is $1\pm\sqrt a \epsilon^{1\over2}$, so that the correction is a half order!
