The answer can be found in <PERSON>'s book Perturbation theory for linear operators. I will use <PERSON>'s notation. In fact I will answer a more general question where you have an operator which depends analytically on a parameter $x$ (your $\epsilon$). Let such operator be $T(x)$ and let
$$ T(x) = \sum_{n=0}^{\infty} x ^n T^{(n)} $$
such that the series converges in a neighborhood of $x=0$. I also call $T=T^{(0)}$. In your case you simply have $T^{(n)}=0$ for $n\ge 2$. We seek the perturbation series of an eigenvalue $\lambda$ of $T$. This means that there exist an eigen-projector of $T$, $P$ such that
$$ TP = \lambda P +D, $$
where $D$ is a nilpotent term that may arise from the <PERSON> decomposition. 768
If $m = \mathrm{dim} P$ is the dimension of the range of $P$, $D^m=0$. Note that for non-degenerate eigenvalue ($m=1$) we have necessarily $D=0$. Define also $Q=1-P$ and the reduced resolvent
$$ S = \lim_{z\to \lambda} = Q (T - z)^{-1} Q $$
lastly let's define
$$ S^{(0)} = -P, \ \ S^{(n)} = S^n, \ \ S^{(-n)} = - D^n, \ \mathrm{for}\ n\ge 1. $$
Let $P(x)$ be the eigenprojector of $T(x)$ analytically connected to $P$. Then one has the following series:
$$ (T(x) - \lambda) P(x) = D + \sum_{n=1}^{\infty} x^n \tilde{T}^{(n)} $$
with
$$ \tilde{T}^{(n)} = - \sum_{p=1}^{\infty} (-1)^p \sum_{\mathcal{A}} S^{(k_1)} T^{(n_1)} S^{(k_2)} \cdots S^{(k_p)} T^{(n_p)} S^{(k_{p+1})}, $$
where $\mathcal{A}$ corresponds to the indices satisfying the following constraint
$$ \mathcal{A} = \left { \sum_{i=1}^p n_i = n ; \sum_{j=1}^{p+1} k_j = p; n_j \ge 1; k_j \ge -m+1 \right }. $$
In the non-degenerate case ($m=1$) this provides the final answer, i.e.
\begin{eqnarray} \lambda(x) &=& \lambda + \sum_{n=1}^{\infty} x^n \lambda^{(n)} \ \lambda^{(n)} &=& \mathrm{Tr} \tilde{T}^{(n)}. \end{eqnarray}
Note that in this case one must have $D=0$. Moreover taking the trace already kills many terms because of the cyclic property of the trace and noting that $SP = PS = 0$.
To make contact with possibly more familiar expressions, note that for a self-adjoint unperturbed operator $T$, the reduced resolvent should look familiar:
$$ S = \sum_{\lambda_j \neq \lambda} \frac{ |j\rangle \langle j|}{ \lambda_j - \lambda}, $$
where I called here $\lambda_j$ and $|j\rangle$ the eigenvalues and eigenvector of $T$. 262