How to prove the following formula, providing a generalization of the Fransén-Robinson constant?
$$\boxed{
\int_0^\infty \frac{t^k}{\Gamma(t)} dt =
\sum_{n=1}^\infty \frac{n^{k+1}}{n!}
+
(2k)!!\sum_{r=0}^k \frac{(-1)^r}{2^{2r}} \binom{k-r}{r} \int_0^\infty \frac{e^{-x}\log^{k-2r}(x)}{[\pi^2+\log^2(x)]^{k-r+1}}dx =: F_k,
}$$
where $k\in\mathbb Z_{\geq 0}$. Here's a run-through for the first few values of $k$ (the sums $S_k$ have been calculated using Dobiński's formula): setting $\varphi_k:=F_k-S_k$,
$$
\begin{array}
FF_0 = &e+\displaystyle\int_0^\infty \frac{e^{-x}}{\pi^2+\log^2(x)}dx = 2.8077702420285\dots \equiv F&\to \quad \varphi_0\approx 0.089489,\\
F_1 = &2e+2\displaystyle\int_0^\infty \frac{e^{-x}\log(x)}{[\pi^2+\log^2(x)]^2}dx = 5.43181977215\dots &\to \quad \varphi_1 \approx -0.004744, \\
F_2 = &5e +2 \displaystyle\int_0^\infty \frac{e^{-x}(\log^2(x)-\pi^2)}{[\pi^2+\log^2(x)]^3}dx= 13.5797413956\dots &\to \quad \varphi_2 \approx -0.011668, \\
F_3 = &15e + 24 \displaystyle\int_0^\infty \frac{e^{-x}(\log^3(x)-\pi^2 \log(x))}{[\pi^2+\log^2(x)]^4}dx = 40.7762149851\dots &\to \quad \varphi_3\approx 0.001988,\ \dots
\end{array}
$$
One could also add
$$
F_{-1} = (e-1)+\frac 1 2 -\frac 1 \pi \int_0^\infty e^{-x} \arctan\left(\frac{\log(x)}\pi \right) dx = 2.2665345077\dots \quad \to \quad \varphi_{-1} \approx 0.548253.
$$
I have arrived at the formula above by observing that the integrand in $\varphi_0$ is such that
$$\frac{e^{-x}}{\pi^2 +\log^2(x)} = \frac {e^{-x}} \pi \mathcal L\{\sin(\pi t)\}(\log x),$$
where $\mathcal L$ is the Laplace transform, and by guessing that the integrand for $k>0$ should similarly involve the Laplace transform of $t^k \sin(\pi t)$. Same for $k=-1$.
My conjecture checks out numerically, but I'd love to know what a rigorous proof of this would look like!
Best Answer
Your argument can be justified and developed. Another one (which is easier for me) is Hankel's formula $$\frac{1}{\Gamma(s)}=\frac{1}{2\pi i}\int_\lambda z^{-s}e^z\,dz$$ (valid for all $s\in\mathbb{C}$), where the contour $\lambda$ encircles the negative real axis (here and below, the principal value of any complex exponentiation is taken). Now, for any $k\in\mathbb{C}$ with $\Re k>-1$, we can write $$\int_0^\infty\frac{t^k\,dt}{\Gamma(t)}=\frac{1}{2\pi i}\int_0^\infty t^k\int_\lambda z^{-t}e^z\,dz\,dt$$ and, if we deform $\lambda$ so that it encircles the disc $|z|\leqslant 1$ (and not only the negative real axis), we get an absolutely convergent double integral, so that the integrations may be interchanged, and we obtain $$\int_0^\infty\frac{t^k\,dt}{\Gamma(t)}=\frac{\Gamma(k+1)}{2\pi i}\int_\lambda\frac{e^z\,dz}{(\log z)^{k+1}}$$ (a side note: the integral on the LHS converges for $\Re k>-2$, so the formula above holds analytically continued, if additionally $k\neq-1$ of course). One may continue this evaluation even for these "general" $k$, by "squeezing" $\lambda$ close to the real axis.
But let's get back to the case of integer $k\geqslant 0$ (see the end for $k=-1$). Then we have $$\int_0^\infty\frac{t^k\,dt}{\Gamma(t)}=k!\ \underbrace{\operatorname*{Res}_{z=1}\frac{e^z}{(\log z)^{k+1}}}_{=A_k}+\frac{k!}{2\pi i}\underbrace{\int_{\lambda'}\frac{e^z}{(\log z)^{k+1}}}_{=J_k},$$ where $\lambda'$ encircles the negative real axis closely. $A_k$ can be computed using the generating function: $$\sum_{k=0}^\infty A_k t^k=\frac{1}{2\pi i}\sum_{k=0}^\infty\int_C\frac{e^z}{\log z}\left(\frac{t}{\log z}\right)^k dz=\frac{1}{2\pi i}\int_C\frac{e^z\,dz}{\log z-t}=\operatorname*{Res}_{z=e^t}\frac{e^z}{\log z-t}=e^{e^t+t}$$ (where the contour $C$ is chosen with $|t/\log z|<1$ uniformly on it).
Hence $k! A_k=eB_{k+1}$ express in terms of the Bell numbers. As for $J_k$, we take the limit of "closely": $$J_k=\int_0^\infty\big((\log x-\pi i)^{-k-1}-(\log x+\pi i)^{-k-1}\big)e^{-x}\,dx,$$ which can be "simplified" using the binomial formula, and giving the result (for $k\geqslant 0$).
For $k=-1$, one may (better) write the integral as $\int_0^\infty\frac{dt}{\Gamma(1+t)}$ and repeat all the steps.