Perturbative solution to $x^3+x-1=0$

The radius of convergence is the largest radius for which there is a holomorphic continuation of your function on a disc centered on $0$.

Here, the root varies holomorphically with $\epsilon$ as long as the roots stay distinct from each other. This happens when the discriminant, $-27-4\epsilon^3$, vanishes.

You can't define the three roots simultaneously as a holomorphic function of $\epsilon$ around any neighborhood of the three bad values of $\epsilon$, because following any small loop around them switches two of them.

As for the branch you're interested in, the real critical point switches the other two roots, so you don't see a problem here. However, you see something happening around the other two critical points.

Therefore, your radius of convergence is $\rho = 3/2^{2/3} \approx 1.889882$.

If you can bound the values of $x$ on the circle of this radius with a $B>0$ (shouldn't be very difficult to show for example that $|x| \le 2$), the Cauchy integral formulas give you an explicit bound $|a_n| \le B \rho^{-n}$ , which allows you to get an explicit exponential bound on the convergence of $x(1)$.

The Lagrange inversion formula applied to the perturbation problem $$ \epsilon x^3+x-1=0 $$ gives the following alternative result. Here $$ f(x)=\epsilon,\qquad\text{ where }\qquad f(x)\equiv\frac{1-x}{x^3}. $$ The unperturbed equation is $f(1)=0$, therefore $$ x(\epsilon)=1+\sum_{n=0}^\infty \frac{g_n}{n!}\epsilon^n, $$ with $$\begin{aligned} g_n &=\lim_{x\to1}\left(\frac{d}{dx}\right)^{n-1}\left(\frac{x-1}{f(x)}\right)^n\\ &=(-1)^n\left(\frac{d}{dx}\right)^{n-1}x^{3n}\big|_{x=1}=\frac{(-1)^n}{2n+1}\frac{(3n)!}{(2n)!}. \end{aligned}$$ Finally, $$ x(\epsilon)=\sum_{n=0}^\infty \frac{(-1)^n}{2n+1}\binom{3n}{n}\epsilon^n. $$ Using the Stirling approximation, for large $n$ the behavior of the $n$th coefficient is $$ \frac{(-1)^n\sqrt3}{2\sqrt{\pi n}(2n+1)}\frac{3^{3n}}{2^{2n}} $$ so that the convergence radius is $$ \frac{2^2}{3^3}. $$ This number is less than $1$ which means that as $\epsilon\to1$, the perturbation series is divergent. But since it is alternating and the absolute value of the coefficients satisfies the Carleman bound $n! 3^n$, the corresponding Padé sequence is ensured to converge to the right answer.

Similar considerations apply to $$ x^3+x-\epsilon=0 $$ where $$ f(x)=\epsilon,\qquad f(x)\equiv x(x^2+1) $$ and $f(0)=0$. Thus, $$ x(\epsilon)=\sum_{n=1}^\infty \frac{g_n}{n!}\epsilon^n $$ with $$\begin{aligned} g_n &=\lim_{x\to 0}\left(\frac{d}{dx}\right)^{n-1}\left(\frac{x}{x(1+x^2)}\right)^n\\ &=\left(\frac{d}{dx}\right)^{n-1}\sum_{l=0}^\infty\binom{n-1+l}{l}(-1)^l x^{2l}\bigg|_{x=0}\\ &=\begin{cases} 0 &\text{if } n-1\text{ is odd}\\ (-1)^k\frac{(3k)!}{k!} &\text{if }n-1=2k\text{ for }k=0,1,2\ldots\ . \end{cases} \end{aligned}$$ To sum up $$ x(\epsilon)=\sum_{k=0}^\infty \frac{(-1)^k}{2k+1}\binom{3k}{k}\epsilon^{2k+1}. $$