There is a slight ambiguity in your question since you do not explain what kind of formula you would like to see: however, since it seems, from the references you cited, that you mean a power series formula, you can build it by using Newton's generalization of the binomial theorem
$$
t\mapsto (1+t)^\alpha = \sum_{k=0}^\infty \binom{\alpha}{k} t^k.\label{1}\tag{1}
$$
By using it, you have that
$$
\left(\sum_{k=0}^\infty a_kx^k\right)^{\alpha}=
\begin{cases}
\displaystyle \sum_{n=0}^\infty \binom{\alpha}{n} \left[\sum_{k=1}^\infty a_kx^k -1\right]^n & a_0=0\\
\\
\displaystyle a_0^\alpha\sum_{n=0}^\infty \binom{\alpha}{n} \left[\sum_{k=1}^\infty a_kx^k\right]^n & a_0\neq0\\
\end{cases}\label{2}\tag{2}
$$
and you could further expand this formula by using the generalization of Cauchy Product formula. However, the computational work involved in the whole procedure seems daunting.
The convergence properties of \eqref{2} can be inferred from the convergence properties of \eqref{1}, listed in the related wikipedia entry, and from the properties of the power series $s(x) = \sum_{k=1}^\infty a_k x^k$: call $D(R)$ the convergence disk (or interval, if $x, a_k\in\Bbb R$ for all $n\in\Bbb N$) of $s(x)$, and $C\big(s^\alpha(x)\big)$ the convergence region of \eqref{2}, i.e the largest open set in $D$ where \eqref{2} converges. Since \eqref{1} converges absolutely for any real $\alpha$ if $|t|<1$ then
$$
C\big(s^\alpha(x)\big) =
\begin{cases}
\big\{x\in D(R) : |s(x)-1|<1\big\} & a_0=0 \\
\\
\big\{x\in D(R) : |s(x)-a_0|<{a_0}\big\} & a_0\neq0 \\
\end{cases}
$$
In a similar fashion, you can analyze the boundary of $C\big(s^\alpha(x)\big)$ by the convergence properties of \eqref{1} on the boundary $|t|=1$ of its convergence disk.
Edit: comments and further notes
Following a comment of the Asker, let's prove the following formula
$$
s^n(x)=\left(\sum_{k=1}^{\infty} a_k x^k\right)^n = \sum_{k=1}^{\infty} \sum_{\substack{0< k_1, \ldots, k_n\le k \\
k_1+\ldots+k_n=k}}a_{k_1} \cdots a_{k_n} x^k \quad\forall n\in\Bbb N,\, n\ge 2\label{3}\tag{3}
$$
Let's proceed by induction
- For $n=2$ equation \eqref{3} is simply the Cauchy Product of the power series $s(x)$ by itself.
- Assuming that \eqref{3} holds true for a general $n\in\Bbb N$ with $n\ge 2$, lets see if it holds for $n+1$. By using Cauchy's product formula,
$$
\begin{split}
\left(\sum_{k=1}^{\infty} a_k x^k\right)^{n+1} & = \left(\sum_{k=1}^{\infty} a_k x^k\right)\cdot\left(\sum_{k=1}^{\infty} a_k x^k\right)^n \\
&= \Bigg( \sum_{k=1}^{\infty} a_k x^k\Bigg)\Bigg(\sum_{k=1}^{\infty} \sum_{\substack{0< k_1, \ldots, k_n\le k \\
k_1+\ldots+k_n=k}}a_{k_1} \cdots a_{k_n} x^k\Bigg)\\
&= \sum_{k=1}^{\infty} \Bigg(\sum_{m=1}^{k} a_m \sum_{\substack{0< k_1, \ldots, k_n\le k-m \\
k_1+\ldots+k_n=k-m}}a_{k_1} \cdots a_{k_n} \Bigg)x^k\\
&= \sum_{k=1}^{\infty} \Bigg(\sum_{\substack{0< k_1, \ldots, k_n, k_{n+1}\le k \\
k_1+\ldots+k_n + k_{n+1}=k}}a_{k_1} \cdots a_{k_n} a_{k_{n+1}}\Bigg)x^k
\end{split}
$$
and the validity of the formula for $n$ implies its validity for $n+1$.
Thus formula \eqref{3} is proven true for every $n\in\Bbb N$ by induction.
A final note. There are several ways of representing the composition of functions $\left(\sum_{k=0}^\infty a_kx^k\right)^{\alpha}$ for a real $\alpha$: for example it is possible to express it by using Faa di Bruno's formula in its various flavors. As stated in the comments, If you know the derivative of the composition of two analytic functions you know the Taylor series of their composition since this is likewise an analytic function. In our specific case we have
$$ s^\alpha(x)=\sum_{k=0}^\infty \frac{1}{k!}\frac{\mathrm{d}^k}{\mathrm{d}x^k}s^\alpha(x)\big|_{x=0} x^k $$
and, as remarked also in the Wikipedia entry on the said formula, you can calculate the Taylor expansion coefficients directly.
What is the "most convenient" way of proceeding? It depends heavily on the structure of the power series $s(x)$ and on on your purpose: in a word it depends on the context of your application.
