In this post the OP proposed the following "Knuth integral" without a rigorous proof $$\int_0^\infty\left(\frac{x^xe^{-x}}{\Gamma(x+1)}-\frac1{\sqrt{2\pi x}}\right)dx=-\frac13$$
I know a proof of corresponding Knuth series by Borwein, which used additional nontrivial techniques other than Lagrange inversion and Lambert W identities, mentioned in the post. However method for that proof is not suitable for this problem.
How can we prove it? Any help will be appreciated.
Here is a solution based on some assumptions. Write
$$ I = - \int_{0}^{\infty} f(x) \, \mathrm{d}x, \qquad f(x) = \frac{x^xe^{-x}}{\Gamma(x+1)}-\frac{1}{\sqrt{2\pi x}}. $$
Adopting the principal branch cut of the complex logarithm, $f$ extends to a holomorphic function on $\mathbb{C}\setminus(-\infty, 0]$. Now we will assume:
Assumption. The line integral of $|f|$ over the circle of radius $R$ vanishes as $R\to\infty$.
This is supported by numerical observations, and I think we may prove this by the Stirling's approximation together with the Euler's reflection formula.1) But I feel exhausted at this moment, so let me leave this part for future updates.
Returning to the computation, we can replace the contour of integration in $I$ by $-x \pm 0^+i$ for $x>0$ by considering upper/lower circular contours. This yields
$$ I = \int_{0}^{\infty} f(-x\pm 0^+i) \, \mathrm{d}x = \int_{0}^{\infty} \left( \frac{e^{\mp i\pi x}x^{-x}e^x}{\Gamma(1-x)} - \frac{\mp i}{\sqrt{2\pi x}} \right) \, \mathrm{d}x. $$
Averaging these two representations and applying the Euler's reflection formula,
$$ I = \int_{0}^{\infty} \frac{\cos(\pi x)x^{-x}e^x}{\Gamma(1-x)} \, \mathrm{d}x = \int_{0}^{\infty} \frac{\sin(2\pi x)\Gamma(x)x^{-x}e^x}{2\pi} \, \mathrm{d}x. $$
Then by applying the integral definition of the gamma function, we get2)
\begin{align*} I &= \int_{0}^{\infty} \frac{\sin(2\pi x)e^x}{2\pi} \left( \int_{0}^{\infty} t^{x-1}e^{-xt} \, \mathrm{d}t \right) \, \mathrm{d}x \\ &= \int_{0}^{\infty} \frac{1}{2\pi t} \left( \int_{0}^{\infty} \sin(2\pi x)(te^{1-t})^x \, \mathrm{d}x \right) \, \mathrm{d}t \\ &= \int_{0}^{\infty} \frac{\mathrm{d}t}{t((2\pi)^2+(1-t+\log t)^2)} \\ &= \int_{-\infty}^{\infty} \frac{\mathrm{d}u}{(2\pi)^2+(1+u-e^u)^2} \tag{$u=\log t$} \end{align*}
To evaluate this integral, let $C_R$ be the contour traversing the square with vertices $\pm R \pm 2\pi i$ in counter-clockwise direction and consider the integral
$$ \int_{C_R} \frac{\mathrm{d}z}{1+z-e^{z}} $$
Using the observation that the equation $1+z-e^{z} = 0$ has a unique solution $z=0$ on the infinite strip $\mathbb{R}\times[-2\pi,2\pi]$, we have
$$ \int_{C_R} \frac{\mathrm{d}z}{1+z-e^{z}} = 2\pi i \, \underset{z=0}{\mathrm{Res}} \left\{ \frac{1}{1+z-e^z} \right\} = \frac{4\pi i}{3}. $$
On the other hand, letting $R\to\infty$ shows that
\begin{align*} \lim_{R\to\infty} \int_{C_R} \frac{\mathrm{d}z}{1+z-e^{z}} = \int_{-\infty}^{\infty} \left( \frac{1}{1+x-e^x - 2\pi i} - \frac{1}{1+x-e^x + 2\pi i} \right) \, \mathrm{d}x = (4\pi i) I. \end{align*}
Therefore we get
$$ I = \frac{1}{3} $$
modulo the assumption.
1) Stirling's approximation shows that $f(z) = \mathcal{O}(|z|^{-3/2})$ as $|z|\to\infty$ along the region $\left|\arg(z)\right|<\pi-\delta$, and I guess that Euler's reflection formula allows to take care of the remaining part by 'reflecting' the gamma function to a more well-behaving region.
2) In fact, Fubini-Tonelli's Theorem is not directly applicable due to the conditionally convergent nature of the double integral. However, this issue can be circumvented by adopting a suitable regularization technique. For instance, inserting the extra factor $e^{-\epsilon x}$ and letting $\epsilon \to 0^+$ will work, since $1-t+\log t \leq 0$.
A solution can be obtained from the formula for $\int_0^\infty\frac{x^x e^{-ax}}{\Gamma(1+x)}\,dx$ given in the linked post (but not proven there), and the asymptotics of Lambert's $\mathrm{W}_{-1}(z)$ as $z\to-1/e$.
In turn, the formula can be obtained from Hankel's integral (also asked for here) $$\frac{1}{\Gamma(s)}=\frac{1}{2\pi i}\int_{-\infty}^{(0^+)}z^{-s}e^z\,dz,$$ where the path of integration encircles the negative real axis (the branch cut of $z^{-s}$).
This is done below (but I prefer to avoid references to "deep properties" of $\mathrm{W}_{-1}$).
For $a>1$ and $x>0$, we put $s=1+x$ and substitute $z=xw$ (the path may be left as is): $$\frac{x^x e^{-ax}}{\Gamma(1+x)}=\frac{1}{2\pi i}\int_{-\infty}^{(0^+)}\left(\frac{e^{w-a}}{w}\right)^x\frac{dw}{w}.$$
Now we deform the contour to have $|e^{w-a}/w|\leqslant c<1$ on it (say, we locate it to the left of $\Re w=b$, and make it enclose $|w|=r$, where $1<r<b<a$). Then integration over $x$ results in an absolutely convergent double integral, and, exchanging the integrations, we obtain ($\lambda$ is the contour described above) $$f(a):=\int_0^\infty\frac{x^x e^{-ax}}{\Gamma(1+x)}\,dx=\frac{1}{2\pi i}\int_\lambda\frac{dw}{w(a-w+\log w)}.$$
The integrand has a pole at $w=w_a$, the solution of $w_a-\log w_a=a$ (here's how Lambert's $\mathrm{W}$ comes in), and the residue is $1/(1-w_a)$; the integral along $|w|=R$ tends to $0$ as $R\to\infty$, giving finally $$f(a)=1/(w_a-1).$$
So, for $a>0$ we have $$\int_0^\infty\left(\frac{x^x e^{-x}}{\Gamma(1+x)}-\frac{1}{\sqrt{2\pi x}}\right)e^{-ax}\,dx=\frac{1}{b}-\frac{1}{\sqrt{2a}},$$ where $b=b(a)$ is the solution of $b-\log(1+b)=a$. The ability to take $a\to 0^+$ under the integral sign is easy to justify (the resulting integral is absolutely convergent, so DCT is applicable).
Thus, the given integral is equal to the limit $$\lim_{b\to 0^+}\left(\frac1b-\frac{1}{\sqrt{2\big(b-\log(1+b)\big)}}\right),$$ which is evaluated easily, and is indeed equal to $-1/3$.
