Prove that $\frac{1}{\phi}<\int_0^\infty \frac{e^{-x}}{\Gamma(x)} dx< \frac{24+\sqrt{2}}{41} $ I'm sure that's a coincidence, but the Laplace transform of $1/\Gamma(x)$ at $s=1$ turns out to be pretty close to the inverse of the Golden ratio: 
$$F(1)=\int_0^\infty \frac{e^{-x}}{\Gamma(x)} dx=0.61985841414477344973$$

Can we prove analytically that: $$\int_0^\infty \frac{e^{-x}}{\Gamma(x)} dx>\frac{1}{\phi}$$

The continued fraction of this number also starts very beautifully:
$$F(1)=\cfrac{1}{1+\cfrac{1}{1+\cfrac{1}{1+\cfrac{1}{1+\cfrac{1}{1+\cfrac{1}{2+\cfrac{1}{2+\cfrac{1}{2+\cfrac{1}{2+\cfrac{1}{2+\dots}}}}}}}}}}=$$
$$=[0; 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 7, 2, 1, 1, 1, 2, 6, 1, 4, 7, 1, 3, 1, 1, 6, 1,\dots]$$
The question has no practical application, but everyone here loves Golden ratio questions, right?


I'm sure there is no closed form for this integral, but if there is some useful transformation, I would like to see it as well.


If we replace every partial quotient in the CF after the first 5 2's by $2$, we obtain a great approximation to the integral:
$$F(1) \approx \frac{1}{41} \left(24+\sqrt{2}\right)$$
So, the more complicated question would be:

Can we prove analytically that: $$\int_0^\infty \frac{e^{-x}}{\Gamma(x)} dx<\frac{1}{41} \left(24+\sqrt{2}\right)$$

The difference is about $4 \cdot 10^{-7}$.
 A: “…if there is some useful transformation, I would like to see it as well.”  Not sure if the following is useful in regard to the problem you pose, however it extends results in the literature, though I’m sure it’s known, (see below) and may possibly help with a bound.
Define
$$F\left( s \right)=\int\limits_{0}^{\infty }{\frac{{{e}^{-sx}}}{\Gamma \left( x \right)}dx}$$
where s is assumed positive and real and then consider the following integral
$$\frac{1}{\Gamma \left( z \right)}=\frac{1}{2\pi i}\int\limits_{-\infty }^{\left( 0+ \right)}{{{e}^{t}}{{t}^{-z}}dt}$$
where the $\left( 0+ \right)$ notation depicts a contour looping around the negative real axis in a positive sense.  Swapping integration leads to
$$F\left( s \right)=\frac{1}{2\pi i}\int\limits_{-\infty }^{\left( 0+ \right)}{\frac{{{e}^{t}}}{s+\log \left( t \right)}dt}$$
Consider now
$$\int\limits_{{{C}_{1}}}^{{}}{\frac{{{e}^{t}}}{s+\log \left( t \right)}dt}$$
where the contour ${{C}_{1}}$ contains a loop integral about the negative real-axis, hooking around the origin connecting to a portion of  a large circle centred at a point $c>{{e}^{-s}}\in {{\mathbb{R}}^{+}}$ which in turn connects to a vertical line in the plane:
$$\int\limits_{{{C}_{1}}}^{{}}{\frac{{{e}^{t}}}{s+\log \left( t \right)}dt}=-\int\limits_{-R}^{\left( 0+ \right)}{\frac{{{e}^{t}}}{s+\log \left( t \right)}dt}+\int\limits_{{{\Gamma }_{R}}}^{{}}{\frac{{{e}^{c+{{\operatorname{Re}}^{it}}}}}{s+\log \left( c+{{\operatorname{Re}}^{it}} \right)}i{{\operatorname{Re}}^{it}}dt}+\int\limits_{\underset{\left( c>{{e}^{-s}} \right)}{\mathop{c-iR}}\,}^{c+iR}{\frac{{{e}^{t}}}{s+\log \left( t \right)}dt}$$
This is a Bromwich contour with a branch cut along the negative-real axis. By construction the contour integral picks up the residue which occurs at $t={{e}^{-s}}$.  Noting that the integral along the semi-circle vanishes as $R\to \infty $
$$\int\limits_{0}^{\infty }{{{e}^{-sx}}dx}\frac{1}{2\pi }\int\limits_{\underset{\left( c>{{e}^{-s}} \right)}{\mathop{-\infty }}\,}^{\infty }{{{e}^{c+it}}{{\left( c+it \right)}^{-x}}dt}=\underset{t={{e}^{-s}}}{\mathop{res}}\,\left( \frac{{{e}^{t}}}{s+\log \left( t \right)} \right)+\frac{1}{2\pi i}\int\limits_{-\infty }^{\left( 0+ \right)}{\frac{{{e}^{t}}}{s+\log \left( t \right)}dt}$$
Note that
$$\underset{t={{e}^{-s}}}{\mathop{res}}\,\left( \frac{{{e}^{t}}}{s+\log \left( t \right)} \right)={{e}^{\frac{1}{{{e}^{s}}}-s}}$$
Then consider $\int\limits_{{{C}_{2}}}^{{}}{{{e}^{c+iz}}{{\left( c+iz \right)}^{-w}}dz}$ where for now $0\le \operatorname{Re}\left( w \right)<1$.  There is a branch point at $z=ic$  and so let ${{C}_{2}}$ have a portion along the real-axis, a semi-circle in the UHP and a branch cut along the imaginary axis incorporating a small circle about the branch point.  As the radius of the semi-circle becomes large and the radius of the circle about the branch point becomes small, the integral over these contours vanishes leaving
$$\int\limits_{-\infty }^{\infty }{{{e}^{c+it}}{{\left( c+it \right)}^{-w}}dt}=2\sin \left( \pi w \right)\int\limits_{0}^{\infty }{{{e}^{-t}}{{t}^{-w}}dt}=2\sin \left( \pi w \right)\Gamma \left( 1-w \right)=\frac{2\pi }{\Gamma \left( w \right)}$$
for $0<\operatorname{Re}\left( w \right)<1$
The integral representation over the real-line is an entire function of w on$\mathbb{C}$and so the result can be extended, for example, to $\operatorname{Re}\left( w \right)\ge 0$.  Therefore
$$\int\limits_{0}^{\infty }{\frac{{{e}^{-sx}}}{\Gamma \left( x \right)}dx}={{e}^{\frac{1}{{{e}^{s}}}-s}}+\frac{1}{2\pi i}\int\limits_{-\infty }^{\left( 0+ \right)}{\frac{{{e}^{t}}}{s+\log \left( t \right)}dt}$$
Breaking up the loop integral, then
$$\int\limits_{0}^{\infty }{\frac{{{e}^{-sx}}}{\Gamma \left( x \right)}dx}={{e}^{\frac{1}{{{e}^{s}}}-s}}+\int\limits_{0}^{\infty }{\frac{{{e}^{-x}}}{{{\pi }^{2}}+{{\left( s+\log \left( x \right) \right)}^{2}}}dx}$$
whereupon making substitutions
$$\int\limits_{0}^{\infty }{\frac{{{e}^{-sx}}}{\Gamma \left( x \right)}dx}={{e}^{\frac{1}{{{e}^{s}}}-s}}+\int\limits_{-\infty }^{\infty }{\frac{{{e}^{x-s-{{e}^{x-s}}}}}{{{\pi }^{2}}+{{x}^{2}}}dx}={{e}^{\frac{1}{{{e}^{s}}}-s}}+\frac{1}{\pi }\int\limits_{-\pi /2}^{\pi /2}{{{e}^{\pi \tan \left( x \right)-s-{{e}^{\pi \tan \left( x \right)-s}}}}dx}$$
This generalises the result in Borwein et. al. “experimentation in mathematics” 2004 pg 287.
Note the following additional perhaps curious/worthless fact:
$${{e}^{\frac{1}{{{e}^{s}}}-s}}=\sum\limits_{n=0}^{\infty }{\frac{{{e}^{-sn}}}{\Gamma \left( n \right)}}$$
and so the whole thing smacks of Euler-Maclaurin.  Going down that path the following pops out
$$\int\limits_{0}^{\infty }{\frac{{{e}^{-st}}\psi \left( t \right)}{\Gamma \left( t \right)}dt}=-0.6198584141460518$$
which to 12 places has the same magnitude as your number $F\left( 1 \right)$.
A: Essentially, this can be proved by dividing the integral interval into $(0,1/100)$, $(1/100,10000)$ and $(10000,\infty)$ and bound the integrand $e^{-x}/ \Gamma(x)$ with simpler functions on each part.
For the lower bound we don't really need to get into the details. There are software that can automatically compute guaranteed numeric upper and lower bounds of integrals. This type of methods are called validated numerics. One such software is arblib. You can access it through SageMath.
The following code proves the lower bound.
low=RBF(CBF.integral(lambda x, _: (e^(-x))/gamma(x), 1/100, 5))
phi=RBF((1+sqrt(5))/2)
low > 1/phi

The output is True.
For the upper bound, 
we use the following inequalities
$$
e^{-x}/\Gamma(x) \le x, \qquad x \in (0,1/100)
$$
and
$$
e^{-x}/\Gamma(x) \le e^{-x}, \qquad x \in (10000,\infty)
$$
So the following code proves the upper bound.
p1=RBF(CBF.integral(lambda x, _: x, 0, 1/100))
p2=RBF(CBF.integral(lambda x, _: (e^(-x))/gamma(x), 1/100,10000))
p3=RBF(integral(e^(-x), x,10000,oo))
up=RBF(1/41*(24+sqrt(2)))
up > p1+p2+p3

The output is again True.

If one want an completely analytic proof, for the lower bound one can use the inequality
$$
f(x):=e^{-x}/\Gamma(x) \ge f_{1,2,5}(x)
$$
where $f_{1,2,5}(x)$ is the Pade approximant of order $(1,2)$ at $x=5$. Since $f_{1,2,5}(x)$ is a rational function, one can do integrate it on $(1/100,5)$ exactly and get an analytic lower bound.
