I would like to calculate the real solution of $$ x^3+x-1=0 $$ by resumming a perturbation series. To this end, I considered $$ x^3+\epsilon x-1=0, $$ $\epsilon$ being a perturbation parameter. The real solution of the unperturbed equation (i.e. for $\epsilon=0$) is $x=1$, so I expanded $x$ in the formal power series $$ x(\epsilon)=1+\sum_{n=1}^\infty a_n \epsilon^n. $$ Substituting this expansion in the above equation and setting to zero the coefficient of $\epsilon^n$ order by order in $n$ gives rise to the following equations $$\begin{aligned} 3a_1+1&=0\\ 3a_2+3a_1^2+a_1&=0\\ 3a_3+6a_1a_2+a_2+a_1^3&=0\\ 3a_4+6a_1a_3+a_3+3a_2^2+3a_1^2a_2&=0 \end{aligned}$$ and so on. Solving recursively for the $a_n$ I got, for for the first few $n$, $$ a_0=1,\ a_1=-\frac{1}{3},\ a_2=0,\ a_3=\frac{1}{3^4},\ a_4=\frac{1}{3^5},\ a_5=0,\\ a_6=-\frac{4}{3^8},\ a_7=-\frac{5}{3^9},\ a_8=0,\ a_9=\frac{77}{3^{13}},\ a_{10}=\frac{104}{3^{14}},\ a_{11}=0,\\ a_{12}=-\frac{595}{3^{17}},\ a_{13}=-\frac{836}{3^{18}}, a_{14}=0,\ a_{15}=\frac{5083}{3^{21}},\ a_{16}=\frac{7315}{3^{22}},\ a_{17}=0. $$ Indeed setting $\epsilon=1$ then $$ x(1)=0,682334\ldots $$ which is a very good approximation to the (explicitly known) real root of the original equation. To justify this approximation, I would like now to calculate the radius of convergence of the perturbation series, and hence I was looking for a closed form of the coefficients $a_n$. However, I couldn't find the general pattern. I also tried looking it up on OEIS but to no avail.
Hints?
EDIT 1: Upon further study, the above coefficients exhibit the following structure. For $n>0$, let us denote $$ a_{n}=-(-1)^{[n/3]}\frac{c_n}{3^{\alpha_n}}. $$ The exponent $\alpha_n$ is given by $$ \alpha_{n}=\sum_{k=0}^\infty\left\lfloor\frac{n}{3^k}\right\rfloor $$ and the nonnegative integer coefficients $c_n$ satsfy $$ c_{3k-1}=0\qquad \text{for }k=1,2,\ldots\ . $$ The first few of them have the following prime number factorizations: ($c_1=c_3=c_4=1)$ $$\begin{aligned} c_6&=2^2\\ c_7&=5\\ c_9&=11 \times 7\\ c_{10}&=13\times 2^3\\ c_{12}&=17\times7\times5\\ c_{13}&=19\times11\times2^{2}\\ c_{15}&=23\times17\times13\\ c_{16}&=19\times11\times7\times5 \end{aligned}$$ In particular, it seems that the same prime number never appears in two consecutive $c_n$. Can we find the general pattern?
EDIT 2: Let us apply the Lagrange inversion theorem, as suggested by Simply Beautiful Art in the comments. Our equation can be rewritten as $$ f(x)=\epsilon,\qquad \text{ where }\qquad f(x)\equiv\frac{1-x^3}{x}, $$ and the unperturbed equation reads $f(1)=0$; then $$\boxed{ x(\epsilon)= 1 + \sum_{n=1}^\infty \frac{g_n}{n!}\epsilon^n } $$ and the coefficients are given by $$\begin{aligned} g_n &= \lim_{x\to 1}\left(\frac{d}{dt}\right)^{n-1}\left(\frac{x-1}{f(x)}\right)^n\\ &=(-1)^n \left(\frac{d}{dx}\right)^{n-1} \left(\frac{x}{1+x+x^2}\right)^n_{x=1}\\ &=(-1)^n \left(\frac{d}{dt}\right)^{n-1} \left(\frac{t+1}{3+3t+t^2}\right)^n_{t=0}. \end{aligned}$$ To calculate these derivatives, let us expand the functions appearing on the right-hand side as follows using the geometric series expansion for small $t$ and Newton's binomial formula: $$\begin{aligned} \left(\frac{t+1}{3+3t+t^2}\right)^n &=\frac{(t+1)^n}{3^n}\frac{1}{\left[1+\left(t+\dfrac{t^2}{3}\right)\right]^n}\\ &=\frac{(t+1)^n}{3^n}\sum_{k=0}^\infty (-1)^k \binom{n-1+k}{k} \left(t+\frac{t^2}{3}\right)^k\\ &=\sum_{k=0}^\infty \sum_{i=0}^k \sum_{j=0}^n \frac{(-1)^k}{3^{n+i}}\binom{n+k-1}{k}\binom{k}{i}\binom{n}{j}t^{k+i+j}. \end{aligned}$$ Now the differentiation is easy: $$ \left(\frac{d}{dt}\right)^{n-1}t^{k+i+j}\bigg|_{t=0}=\frac{(k+i+j)!}{(k+i+j-n+1)!}\delta_{k+i+j,n-1}. $$ We can use the Kronecker delta to get rid of the summation over $j$: we set $j=n-1-i-k$ as long as this value falls within the summation range $$ 0\le n-1-i-k \le n $$ which in turn imposes $i\le n-1-k$ and $k\le n-1$. After a few algebraic simplifications, $$\boxed{ g_n=\sum_{k=0}^{n-1}(-1)^{n-k}(n-1+k)!\sum_{i=0}^{\min(k,n-1-k)}\binom{n}{1+k+i}\frac{1}{3^{n+i}i!(k-i)!}. } $$ Is there a way to compute the radius of convergence of the above series? A rough estimate is given by plugging $k,i\mapsto n/2$ which radius $\rho=2$, but surely a more refined analysis is possible. Numerical evaluations of $1/\lim |g_n/n!|^{1/n}$ suggest that the convergence radius should be strictly less than $2$.
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)$.