Prove the integral evaluates to $\frac{K}{\pi}$
Solution 1:
From a discussion with @Semiclassical on chat I realized I'm going to approach the question totally different and start from the well-known representation of the Catalan's constant, namely $\displaystyle \int_0^{\pi/2} \frac{x}{\sin(x)}=2 K$, and then the right side can be written as
$$\frac{K}{\pi}=\frac{1}{2\pi}\int_0^{\pi/2} \frac{x}{\sin(x)} \ dx=\frac{1}{2}\int_0^{1/2} \frac{\pi x}{\sin(\pi x)} \ dx$$
On the other hand, if we make use of the well-known integral $\displaystyle \int_0^{\infty} \displaystyle \frac{1}{(e^{y}-1)^x}\ dy=\frac{\pi}{\sin(\pi x)}$,
then
\begin{align*}\frac{K}{\pi}&=\frac{1}{2}\int_0^{1/2} x\int_0^{\infty} \displaystyle \frac{1}{(e^{y}-1)^x}\ dy\ dx\\&=\frac{1}{2}\int_0^{\infty}\int_0^{1/2} \displaystyle \frac{x}{(e^{y}-1)^x}\ dx \ dy\\
&=\frac{1}{2}\int_0^{\infty}\left[\frac{1}{\log^2(e^y-1)}-\frac{1}{\sqrt{e^y-1}\log^2(e^y-1)} - \frac{1}{\sqrt{e^y-1}\log((e^y-1)^2)}\right] \ dy\end{align*}
Lastly, by letting $y=x^2$ we're done.
Q.E.D.
Solution 2:
I played a bit around with it Euler style (not worrying to much about establishing convergence; which I think should be fairly simple here if one needed this) and found a way to derive it from scratch.
First make the substitution $t=\log(e^{x^2}-1)$ and then make the second substitution $w=-t$ for the $t<0$ part of the resulting integral to get
$$I = \frac{1}{2}\int_0^\infty \frac{dw}{w^2}\left(1 + e^{-w} - 2e^{-w/2}\right)\frac{1}{1+e^{-w}}$$
We now expand the last part in a Taylor series and introduce (to be able to manipulate each term in the integral independently) a regularization factor $w^s$ into the integrand. This gives us
$$I(s) = \frac{1}{2}\sum_{n=0}^\infty(-1)^n\int_0^\infty dw\left(1 + e^{-w} - 2e^{-w/2}\right)w^{s-2}e^{-nw}$$
where $I(0)$ is the value we are interested in. We can now easily evaluate the integral in terms of the $\Gamma$-function,
$$I(s) = \frac{1}{2}\sum_{n=0}^\infty(-1)^n\left[\frac{1}{n^{s-1}} + \frac{1}{(n+1)^{s-1}} - \frac{2^s}{(2n+1)^{s-1}}\right]\Gamma(s-1)$$
Each term of the sum is well defined in the limit $s\to 0$ and using $\lim\limits_{s\to 0}\Gamma(s-1)s = -1$ together with
$$\lim_{s\to 0} \frac{\frac{1}{n^{s-1}} + \frac{1}{(n+1)^{s-1}} - \frac{2^s}{(2n+1)^{s-1}}}{s} = n\left(2\log(n+1/2)-\log(n+1)-\log(n)\right)\\ + \log(n+1/2)-\log(n+1)$$
gives us
$$I = \frac{1}{2}\sum_{n=0}^\infty(-1)^{n+1}\left[n\left(2\log(n+1/2)-\log(n+1)-\log(n)\right)\\ + \log(n+1/2)-\log(n+1)\right]$$
which after collapsing the logs, turning it into a product and simplifying yields the identity
$$e^{2I-\frac{1}{2}} = \lim_{n\to \infty} \frac{1}{(4n+1)^{2n}}\prod_{k=0}^{n-1}\frac{(4k+3)^{4k+3}}{(4k+1)^{4k+1}}$$
which is identical to the one found here (Eq. 53) with
$$I = \frac{K}{\pi}$$
A proof of the product identity above can be found here.