To derive the identity, it is a straight up application of the Beta function in the form
$$\operatorname{B}(m,n) = 2 \int_0^{\frac{\pi}{2}} \cos^{2m - 1} t \sin^{2n - 1} t \, dt.$$
For the secant integral, we have
\begin{align}
\int_0^{\frac{\pi}{2}} \sec^a t \, dt &= \int_0^{\frac{\pi}{2}} \cos^{-a} t \, dt\\
&= \frac{1}{2} \cdot 2 \int_0^{\frac{\pi}{2}} \cos^{2\left (\frac{1 - a}{2} \right ) - 1} t \sin^{2 \left (\frac{1}{2} \right ) - 1} t \, dt\\
&= \frac{1}{2} \operatorname{B} \left (\frac{1 - a}{2}, \frac{1}{2} \right )\\
&= \frac{1}{2} \frac{\Gamma \left (\frac{1 - a}{2} \right ) \Gamma \left (\frac{1}{2} \right )}{\Gamma \left (\frac{1 - a}{2} + \frac{1}{2} \right )}\\
&= \frac{\sqrt{\pi}}{2 \Gamma \left (1 - \frac{a}{2} \right )} \Gamma \left (\frac{1 - a}{2} \right ),
\end{align}
for $a < 1$, as required. Here the following well-known results of
$$\operatorname{B} (a,b) = \frac{\Gamma (a) \Gamma (b)}{\Gamma (a + b)},$$
together with $\Gamma (\frac{1}{2}) = \sqrt{\pi}$, have been used.
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$.
Best Answer
We only need Euler's formula:
$$\Gamma(1-z) \Gamma(z) = \frac{\pi}{\sin \pi z} \Longrightarrow \Gamma^2\left(\frac{1}{2}\right ) = \pi $$