Finding asymptotic of $ \frac{e^{-n}}{(n-1)!}\int_0^\infty\prod_{k=1}^{n-1}(x+k)e^{-x}dx$ I'm interested in finding the asymptotic at $n\to\infty$ of
$$b_n:= \frac{e^{-n}}{(n-1)!}\int_0^\infty\prod_{k=1}^{n-1}(x+k)\,e^{-x}dx=e^{-n}\int_0^\infty\frac{e^{-x}}{x\,B(n;x)}dx$$
Using a consecutive application of Laplace' method, I managed to get (here)
$$b_n\sim(e-1)^{-n}$$
but this approach is not rigorous, and I cannot find even next asymptotic term, let alone a full asymptotic series.
So, my questions are:

*

*how we can handle beta-function in this (and similar) expressions at $n\to\infty$

*whether we can get asymptotic in a rigorous way ?

 A: First approach. We have
\begin{align*}
b_n & = \frac{{{\rm e}^{ - n} }}{{\Gamma (n)}}\int_0^{ + \infty } {\frac{{\Gamma (x + n)}}{{\Gamma (x + 1)}}{\rm e}^{ - x} {\rm d}x} 
\\ & = \frac{{{\rm e}^{ - n} }}{{\Gamma (n)}}\int_0^{ + \infty } {\frac{1}{{\Gamma (x + 1)}}{\rm e}^{ - x} \left( {\int_0^{ + \infty } {s^{x + n - 1} {\rm e}^{ - s} {\rm d}s} } \right)\!{\rm d}x} 
\\ & = \frac{1}{{\Gamma (n)}}\int_0^{ + \infty } {\frac{1}{{\Gamma (x + 1)}}\left( {\int_0^{ + \infty } {t^{x + n - 1} {\rm e}^{ - {\rm e}t} {\rm d}t} } \right)\!{\rm d}x} 
\\ & = \frac{1}{{\Gamma (n)}}\int_0^{ + \infty } {t^{n - 1} {\rm e}^{ - {\rm e}t} \left( {\int_0^{ + \infty } {\frac{{t^x }}{{\Gamma (x + 1)}}{\rm d}x} } \right)\!{\rm d}t} .
\end{align*}
Employing Ramanujan's formula
$$
\int_0^{ + \infty } {\frac{{t^x }}{{\Gamma (1 + x)}}{\rm d}x}  = {\rm e}^t  - \int_{ - \infty }^{ + \infty } {\frac{{{\rm e}^{ - t{\rm e}^y } }}{{y^2  + \pi ^2 }}{\rm d}y} ,
$$
yields the exact expression
\begin{align*}
b_n & = \frac{1}{{\Gamma (n)}}\int_0^{ + \infty } {t^{n - 1} {\rm e}^{ - ({\rm e} - 1)t} {\rm d}t}  - \frac{1}{{\Gamma (n)}}\int_0^{ + \infty } {t^{n - 1} {\rm e}^{ - {\rm e}t} \int_{ - \infty }^{ + \infty } {\frac{{{\rm e}^{ - t{\rm e}^y } }}{{y^2  + \pi ^2 }}{\rm d}y}\, {\rm d}t} 
\\ & = \frac{1}{{({\rm e} - 1)^n }} - \int_{ - \infty }^{ + \infty } {\frac{1}{{({\rm e} + {\rm e}^y )^n }}\frac{1}{{y^2  + \pi ^2 }}{\rm d}y} .
\end{align*}
Since
$$
\int_{ - \infty }^{ + \infty } {\frac{1}{{({\rm e} + {\rm e}^y )^n }}\frac{1}{{y^2  + \pi ^2 }}{\rm d}y}  \le \frac{1}{{{\rm e}^n }}\int_{ - \infty }^{ + \infty } {\frac{{{\rm d}y}}{{y^2  + \pi ^2 }}}  = \frac{1}{{{\rm e}^n }},
$$
we indeed have
$$
b_n \sim \frac{1}{{({\rm e} - 1)^n }}
$$
as $n\to +\infty$.
Second approach. Changing the order of summation and integration yields
$$
\sum\limits_{n = 1}^\infty  {b_n z^n }  = z\int_0^{ + \infty } {\frac{{{\rm d}x}}{{({\rm e} - z)^{x + 1} }}}  = \frac{z}{{({\rm e} - z)\log ({\rm e} - z)}}
$$
for sufficiently small $z$. Now note that
$$
\frac{z}{{({\rm e} - z)\log ({\rm e} - z)}} = \frac{{\rm e} - 1}{{({\rm e} - 1) - z}} + H(z)
$$
where $H(z)$ is holomorphic in the disc $|z|<\mathrm{e}$. The first term may be expanded as
$$
\frac{{\rm e} - 1}{{({\rm e} - 1) - z}} = \sum\limits_{n = 0}^\infty  {\frac{1}{{({\rm e} - 1)^{n} }}z^n } .
$$
On the other hand, the $n$th Maclaurin series coefficient of $H(z)$ is $\mathcal{O}((\mathrm{e}-\varepsilon)^{-n})$ by the Cauchy–Hadamard theorem for any $\varepsilon>0$ as $n\to+\infty$. Thus
$$
b_n \sim \frac{1}{{({\rm e} - 1)^n }}
$$
as $n\to +\infty$.
A: This is an attempt to generalize the asymptotic. Namely, to consider a bit more general case:
$$b_n(\lambda)=e^{-n}\int_0^\infty\frac{x^\lambda \,e^{-x}}{x\,B(n;x)}dx\,,\,\,\lambda\geqslant0\tag{1}$$
Using the second approach proposed by @Gary, we get
$$\sum\limits_{n = 1}^\infty  {b_n z^n }  = z\int_0^{ + \infty } {\frac{{x^\lambda}}{{(e - z)^{x + 1} }}} dx = \frac{\Gamma(1+\lambda)\,z}{{(e - z)\log^{\lambda+1} (e - z)}}\tag{2}$$
For $\lambda=0, 1, 2,...$ we can obtain the asymptotic in a closed form.
For example, for $\lambda=1$, denoting $\epsilon=e-1-z$ and keeping only the  singular part
$$\frac{z}{{(e - z)\log^2 (e - z)}}=\frac{e-1}{\epsilon^2}-\frac{1}{\epsilon}+H_1(\epsilon)\tag{3}$$
where $H_1(\epsilon)$ is analytical in the disk $\,|\epsilon|<e-1$.
Expanding the singular part of (3) into the series,
$$[z^n]\left(\frac{e-1}{(e-1-z)^2}-\frac{1}{e-1-z}\right)=\frac{1}{(e-1)^{n+1}n!}\big((n+1)!-n!\big)=\frac{n}{(e-1)^{n+1}}\tag{4}$$
It is not clear whether we can find a good approximation for non-integer $\lambda$, though the main asymptotic term for arbitrary $\lambda$ is still $\displaystyle\frac{n^\lambda}{(e-1)^{n+\lambda}}$.
But for $\lambda=k=0, 1, 2, ...$ such closed form does exist. To get it, it is convenient to use the first approach proposed by @Gary.
$$b_n (k)= \frac{{e^{ - n} }}{\Gamma (n)}\int_0^{ + \infty } {\frac{\Gamma (x + n)\,x^ke^{ - x}}{\Gamma (x + 1)} dx}=\frac{\partial^k}{\partial\alpha^k}\,\bigg|_{\alpha=0}\frac{{e^{ - n} }}{\Gamma (n)}\int_0^{ + \infty } {\frac{\Gamma (x + n)\,x^ke^{ - x(1-\alpha)}}{\Gamma (x + 1)} dx}\tag{5}$$
Changing the order of integration and using Ramanujan's formula
$$\frac{{e^{ - n} }}{\Gamma (n)}\int_0^{ + \infty } {\frac{\Gamma (x + n)\,x^k}{\Gamma (x + 1)}e^{ - x(1-\alpha)} dx}=\frac{e^{-n}}{\Gamma(n)}\int_0^\infty e^{-s}s^{n-1}ds\bigg(\int_0^\infty\frac{\Big(se^{-(1-\alpha)}\Big)^x}{\Gamma(1+x)}dx\bigg)$$
$$=\big(e-e^\alpha\big)^{-n}-\int_{-\infty}^\infty\frac{1}{(e+e^{\alpha+y})^n}\frac{dy}{y^2+\pi^2}\tag{6}$$
The second term in (6) is exponentially small vs. the first one.
For example, for $k=1$
$$\frac{\partial}{\partial\alpha}\,\bigg|_{\alpha=0}\int_{-\infty}^\infty\frac{1}{(e+e^{\alpha+y})^n}\frac{dy}{y^2+\pi^2}$$
$$=-n\int_{-\infty}^\infty\frac{1}{(e+e^y)^n}\frac{dy}{y^2+\pi^2}+ne\int_{-\infty}^\infty\frac{1}{(e+e^y)^{n+1}}\frac{dy}{y^2+\pi^2}\sim O\left(\frac{n}{e^n}\right)$$
It mean that the full asymptotics (if we drop exponentially small terms) is given by
$$\boxed{\,\,b_n(k)=e^{-n}\int_0^\infty\frac{x^k \,e^{-x}}{x\,B(n;x)}dx\sim\frac{\partial^k}{\partial\alpha^k}\,\bigg|_{\alpha=0}\big(e-e^\alpha\big)^{-n}\,,\,\,k=0, 1, 2, ...\,\,}$$
where all terms are valid and should be kept.
For several first $k$ we get
$$\qquad\qquad b_n(1)=\frac{n}{(e-1)^{n+1}}\,\,\big(in \,agreement \,with \,(4)\,\big)$$
$$b_n(2)=\frac{n(n+e)}{(e-1)^{n+2}}$$
$$b_n(3)=\frac{n\Big(n^2+3e(n+1)+e^2\Big)}{(e-1)^{n+3}}$$
etc.
The numeric check is in perfect agreement with the approximation. Exactly the same answer can also be obtained by double application of Laplace' method, though in this case we cannot evaluate the error.
