We begin with Stirling's approximation $n!\sim\sqrt{2\pi n}\left(\frac ne\right)^n$; from the fact that $H_n\sim\log(n)+\gamma$, $e^{H_n-\gamma}\sim n$, we can rewrite it $n!\sim\sqrt{2\pi e^{H_n-\gamma}}\left(\frac ne\right)^n$. (This approximation's Laurent series ratio is $1+\frac1{6n}+\ldots$ instead of $1-\frac1{12n}+\ldots$, so this doubles the error and makes it an upper bound, but we need only that it holds.)
So, where $n\rightarrow\infty$, $\sqrt{\frac{2\pi}{e^\gamma}}\sim\frac{n!}{e^{\frac{H_n}2}}\left(\frac en\right)^n=\prod\limits_{k=1}^ne^{1-\frac1{2k}}\frac kn$; however, in order for the partial product to converge to the final value, we can write the $\frac1{n^n}$ as a telescoping product $\prod_{k=1}^n\frac{(k-1)^{k-1}}{k^k}$, which together with the factorial becomes $\prod\limits_{k=1}^n\left(\frac{k-1}k\right)^{k-1}=\frac{n!}{n^n}$, ie.$$\sqrt{\frac{2\pi}{e^\gamma}}\sim\prod\limits_{k=1}^ne^{1-\frac1{2k}}\left(1-\frac1k\right)^{k-1}$$the multiplicands are independent of $n$ (and since $\lim\limits_{k\rightarrow\infty}\left(1+\frac xk\right)^{k+c}=e^x$, they converge to $1$), so the product can be taken to infinity for this limit to become an equality.
Now we take the logarithm (and relabel $k$ to $x$),$$
\frac{\log2\pi-\gamma}2=\sum_{x=1}^\infty1-\frac1{2x}+(x-1)\log\left(1-\frac1x\right)
$$We know $\log\left(1-\frac1x\right)$ (as a complex function of $x$) has poles at $x=0,1$, so its Laurent series $-\sum_{n=1}^\infty\frac1{nx^n}$ converges for $x>1$ (and the $x=1$ case is $0$, since it arises from manipulation that added a $0^0=1$). Multiplying by $x-1$ takes the backward differences of this Laurent series, to $\sum_{n=1}^\infty\frac1n(x^{1-n}-x^{-n})=-1+\sum_{n=1}^\infty\frac1{x^n}\left(\frac1n-\frac1{n+1}\right)=-1+\sum_{n=1}^\infty\frac1{n(n+1)x^n}$, so$$\begin{aligned}
\frac{\log2\pi-\gamma}2&=\sum_{x=1}^\infty1-\frac1{2x}+\sum_{n=1}^\infty\frac1{n(n+1)x^n}\\
&=\sum_{x=1}^\infty\sum_{n=2}^\infty\frac1{n(n+1)x^n}\\
&=\sum_{n=2}^\infty\frac1{n(n+1)}\sum_{x=1}^\infty\frac1{x^n}\\
&=\sum_{n=2}^\infty\frac{\zeta(n)}{n(n+1)}
\end{aligned}$$
Using the polygamma function $\psi_{n-1}(1)=\begin{cases}\gamma&n=1\\(-1)^n(n-1)!\zeta(n)&\mathrm{else}\end{cases}$, this can be written $\log\sqrt{2\pi}=\sum_{n=1}^\infty\frac{\psi_{n-1}(1)(-1)^n}{(n+1)!}$.
It happens that (using Wilf's coeff-extraction notation) $\psi_{n-1}(c)=\left[\frac{x^n}{n!}\right]\log\Gamma(x+c)$, so in this case $\log\sqrt{2\pi}=\sum_{n=1}^\infty(-1)^n\frac1{n+1}\left[x^n\right]\log\Gamma(x+1)=\sum_{n=1}^\infty\frac1{n+1}\left[x^n\right]\log\Gamma(1-x)$. Since turning each $x^n$ into $\frac1{n+1}$ is done by integrating from $0$ to $x$,$$\log\sqrt{2\pi}=\int_0^1\log\Gamma(1-x)\ dx$$
To account for $α>0$, we instead transform Stirling's approx. into$$\begin{aligned}
\sqrt{\frac{2\pi}{e^\gamma}}\left(\frac{αe^\gamma}e\right)^α&\sim n!e^{n+(α-\frac12)H_n}\frac{α^α}{(α+n)^{α+n}}\\
&=\prod_{k=1}^ne^{1+\frac{α-\frac12}k}\frac k{α+k-1}\left(1-\frac1{α+k}\right)^{α+k}
\end{aligned}$$giving$$\begin{aligned}
\log\left(\sqrt{\frac{2\pi}{e^\gamma}}\left(\frac{αe^\gamma}e\right)^α\right)&=\sum_{k=1}^\infty 1+\frac{α-\frac12}k+\log\left(\frac k{α+k-1}\right)+(\alpha+k)\log\left(1-\frac1{α+k}\right)\\
&=\sum_{k=1}^\infty 1+\frac{α-\frac12}k+\left(\sum_{i=1}^\infty\frac{\left(\frac{1-\alpha}x\right)^i}i\right)+\left(\fracαx+1\right)\left(\sum_{i=0}^\infty\frac{α^i}{(-x)^{i+1}}\sum_{j=0}^i{i\choose j}\frac{α^{-j}}{j+1}\right)\\
&=\sum_{x=1}^\infty\frac{α-\frac12}x+\left(\sum_{i=1}^\infty\frac{\left(\frac{1-\alpha}x\right)^i}i\right)+\left(\sum_{i=1}^\infty\left(-\frac αx\right)^i\sum_{j=1}^i{i-1\choose j-1}\frac1{(j+1)(-a)^j}\right)\\
&=\sum_{x=1}^\infty\sum_{i=3}^\infty\frac1{i(i-1)x^{i-1}}\sum_{j=0}^{i-1}{i\choose j}(-\alpha)^j\\
&=\sum_{i=3}^\infty\left((1-\alpha)^i-(-\alpha)^i\right)\left(\frac{\zeta(i-1)}{i(i-1)}=\frac{\psi_{i-2}(1)}{i!}=\frac{\left[x^{i-1}\right]\log\Gamma(1+x)}i\right)\\
&=\int_{\alpha-1}^\alpha\log\Gamma(1+x)\ dx+\int_{\alpha-1}^\alpha\gamma x\ dx\\
&=\int_0^1\log\Gamma(\alpha+x)\ dx+\left(\alpha-\frac12\right)\gamma
\end{aligned}$$
ie. $\log\left(\sqrt{2\pi}\left(\fracαe\right)^α\right)=\int_0^1\log\Gamma(\alpha+x)\ dx$ as requested.
Since each step is two-way, performing this derivation in reverse recovers (an also-valid asymptotic that can be simplified to) Stirling's approximation from both the $\alpha=0$ case and the general Raabe's formula.
Note that this answer (on a duplicate of this question) is much simpler and shows it directly via Gauss's multiplication theorem (which is also reversible, by using the gamma function's log-convexity), and simple arguments via limits obtain Stirling's approximation and Gauss's theorem from each other.
