Is it possible to derive an asymptotic expression as $n \to \infty$ for the following sum:
$$S_n(a,x) = \sum_{k=1}^n \binom{n}{k} \frac{k!^2 x^k}{k \cdot (a)_k} $$
$(a)_k$ is the rising factorial (Pochhammer symbol).
Euler-Maclaurin seems way too complicated here, and I'm not sure what else to try.
Einar Rødland showed that $$ S_n (a,x) \approx n!B(a,n)x^n {}_1F_1\! \left( { - n - a + 1, - n + 1;\tfrac{1}{x}} \right). $$ Here $B$ is the beta function and ${}_1 F_1$ is one of Kummer's hypergeometric functions. Employing Kummer's transformation ${}_1F_1 (- n - a+1;- n+1;\tfrac{1}{x}) = e^{1/x} {}_1F_1 (a;1 - n; - \tfrac{1}{x})$ and truncating the series for ${}_1F_1$ after $n$ terms, we obtain the approximation $$ S_n (a,x) \approx \frac{{n!}}{{\Gamma (n + a)}}x^n e^{1/x} \sum\limits_{k = 0}^{n - 1} {\frac{{\Gamma (k + a)}}{{k!x^k }}\Gamma (n - k)} . $$ Note that, in general, the terms in the sum are decreasing in magnitude as $k$ is increasing. Numerically this formula seems to perform very well even when $|a/x|$ is comparable with or larger than $n$. Indeed, the ratio of the last and the first term is $$ \frac{{\Gamma (n + a - 1)}}{{\Gamma ^2 (n)\Gamma(a)x^{n - 1} }}, $$ which is small if $x$ is not close to $0$ even if $a=\mathcal{O}(n)$. For large $n$, one may truncate the series after any finite number of terms to obtain simpler approximations. In particular, to leading order, $$ S_n (a,x) \approx n!B(a,n)x^n e^{1/x} . $$
I followed the same suggestion as @Maxim of summing over $r=n-k$. The main reason is that for large $n$, ie $n\gg a, 1/|x|$, the largest terms are with $k$ close to $n$: having them for small $r$ is more convenient.
This gave me a somewhat different looking expression: $$ S_n(a,x) = \frac{n!(n-1)!(a-1)!x^n}{(n+a-1)!}\cdot \sum_{r=0}^{n-1}\frac{(n+a-1)^{\underline{r}}}{(n-1)^{\underline{r}}r!x^r} = \frac{(n-1)!x^n}{\binom{n+a-1}{a-1}}\cdot \sum_{r=0}^{n-1}\frac{\binom{n+a-1}{r}}{\binom{n-1}{r}r!x^r} $$ where $m^{\underline{r}}=m!/(m-r)!$ is the falling factorial. The sum actually corresponds to a generalised hypergeometric series $$ S_n(a,x) \sim \frac{(n-1)!x^n}{\binom{n+a-1}{a-1}}\cdot {}_1F_1\left( \genfrac..{0pt}{}{-(n+a-1)}{-(n-1)} ;\frac{1}{x}\right) $$ although this does not exist when $n,a$ are positive integers as the denominator becomes zero for $r\ge n$ so the GHS is not naturally truncated to $r<n$.
I'm not sure what type of arguments are of interest. For large $n$, we may approximate $(n+a-1)^{\underline{r}}/(n-1)^{\underline{r}} \approx [(n+a-1)/(n-1)]^r$, which yields $$ S_n(a,x) \approx \frac{(n-1)!x^n}{\binom{n+a-1}{a-1}} \cdot \exp\left(\frac{n+a-1}{(n-1)x}\right). $$ A few test computations where $n\gg a,1/|x|$ looked promising. Maxim's approximation seems to require $n\gg a/|x|$, which this approximation is somewhat less dependent on. Of course, both fail close to $x=0$.
Another improvement in terms of efficient approximation is to use a power series expansion of the logarithm of GHS, since these often converge much faster: eg a low order approximation of $\ln F(-(n+a-1);-(n-1);u)$ may outperform the above approximation which is just to the first order.
To assess the terms of the original sum more easily, we may rewrite it $$ S_n(a,x) = \sum_{k=1}^n \frac{n^\underline{k}(k-1)!x^k}{a^\overline{k}} = \frac{nx}{a}\cdot\sum_{j=0}^{n-1} \frac{(n-1)^\underline{j\,}j!x^j}{(a+1)^\overline{j\,}}. $$ using $p^\overline{q}$ and $p^\underline{q}$ to represent rising and falling factorials: I prefer these over the Pochhammer symbols which may be used to represent either.
In general, the terms (absolute value in case $x$ is negative) may be divided into three regions:
As $x$ approaches 0, the rising mid-range will shrink until, for $x$ close to 0, the terms will be falling throughout. The first term will then be the dominant one, while adding a few more terms may improve the approximation.
As $|x|$ increases, the mid-range rising region expands, while first the low-end falling region and later the high-end falling region shrink. The high-end terms will be dominant and similar to the terms of the exponential series as used in the above approximations.
Difficulties are caused by two issues. First, the terms may have two local maxima: one at $j=0$ ($k=1$) and another between the mid- and high-range regions. Either of these may be the dominant one.
The second issue is that, when the higher term is dominant, the expomential appoximation works better when the dominant term is close to the upper end: ie for $k$ close to $n$, or low $r=n-k$ after. In this case, my original idea for analysing the sum might work better: to identify which $k$ gave the dominant term, and then build the approximation around that.
Note that the ratio between consecutive terms is $(n-k)kx/(a+k)$, so it is fairly easy to check where the terms are rising or falling. Also, solving in $k$ for ratio close to $1$ will help identify the maximum and then approximate nearby terms relative to this.
Not an answer.
$$S_n(a,x) = \sum_{k=1}^n \binom{n}{k} \frac{(k!)^2 }{k \, (a)_k}x^k$$ $$\frac{\partial S_n(a,x)}{\partial x} = \sum_{k=1}^n \binom{n}{k} \frac{(k!)^2 }{ \, (a)_k}x^{k-1}$$ $$\frac{\partial S_n(1,x)}{\partial x}=\frac{e^{\frac{1}{x}}\, x^n\, \Gamma \left(n+1,\frac{1}{x}\right)-1}{x} \quad \implies \quad S_n(1,x)=\color{red}{\large ???}$$
Not an answer (so I still hope for someone to post a better one and get the bounty), but I need somewhere to put my notes on this.
Using @Maxim 's suggestion in the comments, we change the summation index.
Denote $k=n-l$, then
$$S_n(a,x) = n! x^n \sum_{l=0}^{n-1} \frac{(n-l-1)!}{l! ~ (a)_{n-l} ~ x^l}$$
$$(a)_{n-l}=\frac{\Gamma(a+n-l)}{\Gamma(a)}$$
$$S_n(a,x) = \Gamma(a) n! x^n \sum_{l=0}^{n-1} \frac{\Gamma(n-l)}{l! ~ \Gamma(n-l+a) ~ x^l}$$
$$S_n(a,x) = n! x^n \sum_{l=0}^{n-1} \frac{B(n-l,a)}{l! x^l}$$
$$B(n-l,a)=\int_0^1 t^{n-l-1}(1-t)^{a-1}dt$$
However, Beta function is not very useful here. Instead, we could use Stirling approximation for Gamma function, which up to 2 orders, reads:
$$\Gamma(n+1)\asymp {\sqrt {2\pi n}}\left({\frac {n}{e}}\right)^{n}\left(1+{\frac {1}{12n}}+{\frac {1}{288n^{2}}}+O(1/n^3) \right)$$
Thus we have for the 0 order approximation, assuming $l$ fixed:
$$\frac{\Gamma(n-l)}{\Gamma(n-l+a)} \asymp e^a \left(\frac{n-l-1}{n-l-1+a}\right)^{n-l-1/2} \frac{1}{(n-l-1+a)^a} $$
So, we have asymptotic series:
$$S_n(a,x) \asymp \Gamma(a) n^{-a} n! x^n \sum_{l=0}^{m} \frac{1}{l!x^l}$$
Where $m$ is some number, $m\ll n$.
However, I checked numerically for some values, and it's a very poor approximation.
We may need to use more orders of Stirling approximation. I'll return to this problem as soon as I'm able.