Given the perturbation algebraic problem $z^2(1-(\varepsilon z) + (\varepsilon z)^2) +1 = 0 $ solutions when $\varepsilon \rightarrow 0$ , I wish to evaluate the first order terms of the 4 roots.
This singular perturbation problem is caused by the fact that there are 2 solutions $r_1 ,r_2 $ such that $|r_1| , |r_1| \approx 1$ and the other $2$ are of order $|r_3|, |r_4| \approx \varepsilon^{-1}$. So I wish to find $p\in \mathbb{R}$ values such that $z(\varepsilon) = \varepsilon^{-p} \sum_1^\infty a_i\varepsilon^i$ to approximate all 4 roots, in particular to find $a_{0_1}, a_{0_2},a_{0_3}, a_{0_4}$ the first term of the approximation.
Rearranging the equation : $\varepsilon^2z^4 - \varepsilon z^3 + z^2 +1 = 0$ , and assigning $z=\varepsilon^{-p}\sum_0^\infty a_i\varepsilon^i $ we get: $\varepsilon^{2-4p}(\sum_1^\infty a_i\varepsilon^i)^4 - \varepsilon^{1-3p}(\sum_1^\infty a_i\varepsilon^i)^3 + (\sum_1^\infty a_i\varepsilon^i)^2 +1 =0$
to for $p\le0$ taking the limit $\varepsilon \rightarrow0$ we are left with $a_0^2 +1=0\Rightarrow a_0\in\{i,-i\}$. So we found first order approximation for $r_1,r_2$ .
Now in order to get different limits we need $1-3p < 0$ or $2-4p <0$ so $\varepsilon^{1-3p}$ or $\varepsilon^{2-4p}$ tend to infinity when $\varepsilon \rightarrow 0$ , so we we are left to examing $p\ge \frac{1}{3}$.
When $p = \frac{1}{3}$ and $\varepsilon \rightarrow 0$ the equation is: $0 -a_0^3+a_0^2+1=0$ this result in $3$ approximation such that $|a_0| \approx 1$ and we are searching for $|a_0| \approx \varepsilon^{-1}$ i.e "large" roots.
when $\frac{1}{3} < p < \frac{1}{2}$ and $\varepsilon \rightarrow 0$ the equation is : $lim_{\varepsilon \rightarrow 0} -\varepsilon^{1-3p}a_0^3+a_0^2 +1=0$ let us neglect $1$ because when $\varepsilon \rightarrow 0$ $\varepsilon^{1-3p}\rightarrow\infty$ so we are left with $a_0^2(-\varepsilon^{1-3p}a_0+1)=0$ so $a_0 = \varepsilon^{-1+3p}$ is a first order term for approximating one of the big roots. Yet non of those $p$ values gives us two values of $a_0$. So I tried to search for $p\ge 1/2$ yet, the same thing is going to happen, let us show it quickly:
taking $p\ge 0$ and $\varepsilon \rightarrow 0$ we get: $\varepsilon^{2-4p}a_0^4 -\varepsilon^{1-3p}a_0^3 + a_0^2 + 1 = 0$ , again neglecting $1$ and dividing by $a_0^2$ er are left with $\varepsilon^{2-4p}a_0^2 -\varepsilon^{1-3p}a_0 + 1 = 0$ let us again neglect $1$ from the same reasons: $\varepsilon^{2-4p}a_0 -\varepsilon^{1-3p}=0$ so $a_0 = \varepsilon^{1-3p -2+4p} = \varepsilon^{-1+p} $. Which is again just $1$ root.
I wish to get a specific $p$ which gives two possible $a_0$ values as the first term for the large roots. One mistake I might have done is to neglect the possibility for $a_0 = 0$ (when dividing by $a_0$) , however in this case that $a_0=0$, it seems that for all $p\ge \frac{1}{3}$ $z(\varepsilon) = \varepsilon^{-p} \sum_0^\infty a_i\varepsilon^i$ gives an approximation for both large roots, in this case I don't understand why there is no unique value for $p$?
Edit: after receiving the answer, I figured out that I made a mistake (can be seen above) when assigning $z =\varepsilon^{-p}\sum_0^\infty a_i\varepsilon^i$ to the equation, I forgot to multiply by $\varepsilon^{-2}$ the $z^2$ term. However I wish to show why $p=-1$ is the only $p$ which gives an approximation $\varepsilon^{-p}\sum_0^\infty a_i\varepsilon^i$ for the large roots. So assigning the approximation series when $p>0$ to the equation we (this time correctly) get: $\varepsilon^{2-4p}(\sum_1^\infty a_i\varepsilon^i)^4 - \varepsilon^{1-3p}(\sum_1^\infty a_i\varepsilon^i)^3 + \varepsilon^{-2p}(\sum_1^\infty a_i\varepsilon^i)^2 +1 =0$ let us write $p=1+q$ for $q\in(-1,\infty)$
So: $\varepsilon^{-2-4q}(\sum_1^\infty a_i\varepsilon^i)^4 - \varepsilon^{-2-3q}(\sum_1^\infty a_i\varepsilon^i)^3 + \varepsilon^{-2q-2}(\sum_1^\infty a_i\varepsilon^i)^2 +1 =0$ we are intrested only in the smallest powers of $\varepsilon$ when taking $\varepsilon \rightarrow 0$ so we may neglect any coefficient which is not $a_0$ to get: $\varepsilon^{-2-4q}a_0^4 - \varepsilon^{-2-3q}a_0^3 + \varepsilon^{-2 -2q}a_0^2 =0 \Rightarrow \varepsilon^{-2q}a_0^2 - \varepsilon^{-q}a_0 + 1 =0 $ in case $q\neq 0$ $\varepsilon^{-2q}$ and $\varepsilon^{-2q}$ tend to infinity so we neglect $1$ and left only with one $a_0$ value. So $q$ must be $0$ in order to get quadratic equation which gives $a_{0_{1,2}} = \frac{1\pm\sqrt{3}i}{2}$.
The first order approximations of the four roots of this quartic equation can be found by allowing the equation to degenerate in two different ways. You already figured out the first one: simply setting $\epsilon=0$ yields the quadratic equation $z^2+1=0$, whose two roots $r_1,r_2=\pm i$ are the first order approximations you already found.
However, there is another way to degenerate the quartic that you didn't use. Namely, change variables to $x=\epsilon z$ (this scaling makes sense if we are searching for roots of order $\epsilon^{-1}$) yielding the equation $$ (x/\epsilon)^2(1-x+x^2)+1=0, $$ or equivalently $$ 1-x+x^2+\epsilon^2 x^{-2}=0 $$ (the division is valid since $x=0$ is not a root). Setting $\epsilon$ to $0$ yields the second degeneration $$ 1-x+x^2=0,\qquad x=\frac{1\pm i\sqrt{3}}{2} $$ by the quadratic formula. Since $z=x/\epsilon$ we obtain the first order expressions for the other two roots, $$ r_3,r_4=\frac{1\pm i\sqrt{3}}{2\epsilon}. $$ It is interesting to note that $r_1,r_2$ are a conjugate pair of $4^{th}$ roots of unity, whereas $r_3,r_4$ are a conjugate pair of $6^{th}$ roots of unity (rescaled by $\epsilon^{-1}$).