root perturbation of a polynomial and the method of dominant balance

Consider the following polynomial equation $f(z)=z^{2m}+a(1+z)^m=0$, with $a\to 0$. A simple argument using Rouche's theorem shows that $z_k = (1+o(1))\zeta_k a^{1/2m}$, where $\{\zeta_k\}$ are roots of $\zeta^{2m}+1=0$.

I would like to obtain higher-order approximation of $z$. In particular, I wish to show that $z_k = \zeta_k(1) a^{1/2m} + \zeta_k(2) a^{1/m} + \cdots + \zeta_k(L)a^{L/2m} + o(a^{L/2m})$, with coefficients $\zeta_k(1),\cdots,\zeta_k(L)$ only depending on $2m$ but not $a$.

I think by using the so-called method of dominant balance I can derive the values of the coefficients. This document

provides an easy introduction to this method. Basically, by using $L$-th order approximation and balancing the coefficients carefully, one can obtain $\tilde z_k=\xi_k(1)a^{1/2m}+\cdots+\zeta_k(L)a^{L/2m}$ that satisfy $f(\tilde\zeta_k) = O(a^{1+\frac{L}{2m}})$.

My question is: how to prove that such a power seris expansion exists in the first place? I searched for a long time but it seems most solution assumes such an expansion at the first place and then proceeds to matching the terms. Note also that $|f(\tilde z_k)|$ being small does not imply that $\tilde z_k$ is close to a root of $f(z)=0$.

More generally, let $P_a(z)$ be a polynomial of degree $m$, with its coefficients analytical functions of a parameter $a$ that goes to 0. Under what conditions can we conclude that all roots of $P_a(z)=0$ admit an analytical expansion of $z_k = \sum_{j=0}^{\infty}{\zeta_k(j) a^{j/m}}$, in a small neighborhood of roots of $P_0(z)=0$?

If $z = \zeta s v$ with $\zeta$ a $2m$-th root of $-1$ and $s = a^{1/(2m)}$, the equation becomes $$1 + \zeta s v = v^2$$ and then $$v = \frac{\zeta s}{2} - \sum_{k=0}^\infty {2k \choose k} \frac{(-1)^k}{2k-1} \left(\frac{\zeta s}{4}\right)^{2k}$$