** About the non-convergence
Non-convergence of the formal powerseries does not bother
Several times it was mentioned in the answers that the formal powerseries of the half iterate does not converge. **
There That is true, however there is an elaborated theory about the fractional iteration of analytic functions with a fixed point $z_0$. z_0$ which gives more far reaching answers.
Here we have the case of a parabolic fixpoint, i.e. $f'(z_0)=1$. These functions have mostly no fractional iterates analytic at the fixpoint.
But They have unique fractional iterates to the sides of the fixpoint, i.e. there are several domains bounded by/around the fixpoint which have the formal powerseries as asymptotic powerseries.
The arrangement of these domains is called Leau-Fatou flower (See the online book of Milnor [3] for details). The petals are alternating attractive and repellent when following the circle around the fixpoint. The number of these domains/petals is determined by number $m$ of zeros after the coefficient 1 in the powerseries development of $f$ at $z_0$. The number of domains is $2(m+1)$.
In our case the fixpoint is 0 and the development is $e^z-1=z+\frac{z^2}{2}+\dots$, so $m=0$ and the number of petals is 2. One petal (the repelling) is on the positive axis and one petal (the attracting) is on the negative axis. On these two petals (which overlap in the complex plane) are the two (different, not being analytic continuations) solutions defined, that have the formal powerseries as asymptotic powerseries.
There are several (general) formulas possible to numerically compute these two solutions.
The classic formula of Lévy for the Abel function (with $\alpha_u(u)=0$) is too slow for computations: $$\alpha_u(z) =\lim_{n\to\infty}\frac{f^{[n]}(z) - f^{[n]}(u)}{f^{[n+1]}(u)-f^{[n]}(u)} $$
The Newton formula for the regular fractional iteration is also too slow: $$f^{[t]}(z) = \sum_{n=0}^\infty \binom{t}{n} \sum_{m=0}^n \binom{n}{m} (-1)^{n-m} f^{[m]}(z)$$
But the following formulas for the Abel function (adapted to $f(x)=e^x-1$) are quickly converging: $$\alpha_1(z) = \lim_{n\to\infty} \frac{1}{3}\log(-f^{[n]}(z)) - \frac{2}{f^{[n]}(z)} - n, \quad z<0$$ $$\alpha_2(z) = \lim_{n\to\infty} \frac{1}{3}\log(f^{[-n]}(z)) - \frac{2}{f^{[-n]}(z)} + n, \quad z>0$$
You get the half iterate from the Abel function by $f^{[1/2]}(z)=\alpha^{-1}(1/2+\alpha(z))$ (independent on any additive constant of the Abel function). The non-Lévy formulas are probably first discovered by Écalle in his thesis [2] which deals completely with the parabolic case $f'(z_0)=1$.
[1] Kuczma, M., Choczewski, B., & Ger, R. (1990). Iterative functional equations. Encyclopedia of Mathematics and Its Applications, 32. Cambridge University Press.
[2] Écalle, J. (1974). Théorie des invariants holomorphes. Publications math'ematiques d'Orsay, 67-74 09. Orsay: Univ. Paris-XI.
[3] Milnor, J. (2006). Dynamics in one complex variable. 3rd ed. Princeton Annals in Mathematics 160. Princeton, NJ: Princeton University Press. viii, 304 p.

