The sum in question can be written (using Latin letters instead of Greek) as $$ S_n(a,b):= \sum_{k=0}^n {n \atopwithdelims \{ \} k} k! \binom{na}{k} b^k $$ where the round brackets denotes an ordinary binomial. It will be shown as $n \to \infty$ $$ (1)\quad S_n(a,b) \sim \frac{n!}{2}\,\exp{\Big(n\big( h(u_0) + \frac{h''(u_0)}{2}\,u_0^2\big)\Big)} \, \text{erfc}\big(\sqrt{\frac{n}{2} h''(u_0)} \,\,u_0 \big) $$ where $$h(u) = \frac{\log(1+b(e^u-1))^a}{u}, \quad u_0=\frac{1}{a} + \text{W}\big(\frac{1-b}{a\,b}\,e^{-1/a}\big) $$ and $W()$ is the Lambert W function. This follows from an asymptotic analysis of an alternative form,
$$(2) \quad S_n(a,b)=n!\,[u^n] (1+b\,(e^u-1))^{n \,a}$$
where $[u^n]$ is the 'coefficient of' operator. Note that when $b=1,$ the sum can be solved exactly, and gives the form the OP suspects. (This summation is known.) To derive the alternative form, note that sum is the Hadamard product of
$$ \omega_n(z) = \sum_{k=0}^\infty {n \atopwithdelims \{ \} k} k! z^k \quad \text{& } \quad (1+z)^{na}=\sum_{k=0}^\infty \binom{na}{k} z^k $$ The expression for $\omega_n$ can be written in terrms of Eulerian polynomials, but I prefer the polylogarithm at negative argument, $$ \omega_n(t) = \frac{\text{Li}_{-n}(t/(1+t))}{1+t} $$ Taking the Hadamard product,
$$S_n(a,b) = \frac{1}{2 \pi \,i} \oint \frac{dt}{t(1+t)} \text{Li}_{-n}(\frac{t}{1+t}) \big(1+\frac{b}{t}\big)^{na} \, dt$$ where the integration path surrounds the origin. Use $\text{Li}_{-n}(t/(1+t))=(-1)^{n+1}\text{Li}_{-n}(1+1/t).$ Let $t \to 1/t$ in the integral and check residue at $\infty.$ Then $$S_n(a,b) = \frac{(-1)^n}{2\pi i}\oint \frac{ \text{Li}_{-n}(1+t)}{1+t} \big(1+b\,t\big)^{n\,a} \, dt$$ Now use the known generating function $$\sum_{n=0}^\infty \frac{u^n}{n!} (-1)^n \frac{ \text{Li}_{-n}(1+t)}{1+t} = \frac{1}{e^u-1-t}$$ Putting this equation in the penultimate and using Cauchy's theorem results in equation (2).
For the asymptotic analysis, write (2) as a countour integral $$S_n(a,b)=\frac{n!}{2 \pi \,i} \oint \Big( \frac{(1+b(e^u-1))^a}{u} \Big)^n \frac{du}{u} $$ Classic saddle point analysis tells us to write this as $$S_n(a,b)=\frac{n!}{2 \pi \,i} \oint \exp{\big(n h(u)\big)} \frac{du}{u} $$ where $h$ has been given in (1), then expand $h$ in a power series about $u_0,$ which has also been found explicitly and given in (1). Run the contour through $u_0$ and open the loop so that it becomes a vertical line. Then, so long as certain conditions are satisfied (I have not checked them, other than $h''(u_0)>0$, a necessity) then $$ S_n(a,b) \sim \exp{(n\ h(u_0))} \frac{n!}{2 \pi} \int_{-\infty}^\infty \frac{dt}{u_0+i\,t} \exp{\big(-\frac{n}{2}\ h''(u_0) t^2 \big)} $$ The integral is solvable in terms of the complementary error function. For $n$ sufficiently large, it will simplify in the asymptotic limit to an exponential, as long as $h''(u_0)$ isn't too small. A potential problem with this analysis is that if $h''(u_0)$ does become small, and much smaller than $h^{(3)}(u_0)$ then a quadratic expansion of $h$ is not sufficient.
For examples of the accuracy of the approximation, a table is shown. $$ \begin{array}{c|lcr} n & a & b & \text{true} & \text{asym} & \text{% err}\\ \\ \hline 20 & 1/2 & 1/2 & 2.7717\ 10^{17} & 2.6867\ 10^{17} & 3.07\\ {} & 1/2 & 3/2 & 3.7421\ 10^{21} &3.5745\ 10^{21} & 4.48 \\ {} & 3/2 & 1/2 & 2.1498\ 10^{25} & 2.0862\ 10^{25} &2.96\\ {} & 3/2 & 3/2 & 1.6663\ 10^{23} & 1.5884\ 10^{23} & 4.68\\ 200 & 1/2 & 1/2 & 6.634\ 10^{374} & 6.617\ 10^{374} & 0.34\\ {} & 1/2 & 3/2 & 3.336\ 10^{415} & 3.319\ 10^{415} & 0.51\\ {} & 3/2 & 1/2 & 6.158\ 10^{453} & 6.138\ 10^{453} & 0.33\\ {} & 3/2 & 3/2 & 8.720\ 10^{521} & 8.673\ 10^{521}& 0.53 \end{array} $$
It's really surprising that the approximation for a problem with 2 parameters is so easily characterized.
