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.