The Gudermannian function and Apéry's constant

You already asked ten similar questions. Expand $\frac{1}{\cosh y}$ in powers of $e^y$, integrate term by term, expand $\frac{1}{1+e^{2y}}$ in powers of $e^y$, multiply the two series, and use that $\int_0^\infty y e^{-ny}\mathrm dy = \int_0^\infty (x/n) e^{-n(x/n)}\mathrm d(x/n)=n^{-2} \int_0^\infty x e^{-x}\mathrm dx$ $ = n^{-2} (-xe^{-x}|_0^\infty +\int_0^\infty e^{-x}\mathrm dx) = n^{-2}$

I would like to present another approach directly utilizing the definition $\operatorname{gd}(x)=2\arctan(e^x)-\frac\pi2$, which is valid for real $x$, as Brevan Ellefsen suggested within the comments.

First of all note that the integrand is an even function from which we can conclude that it is not of matter whether we integrate over the interval $[-\infty,0]$ or $[0,\infty]$ so will use the latter one for the sake of simplicity. Using the well-known relation $\arctan(x)+\arctan\left(\frac1x\right)=\frac\pi2$ we can simplify the given integral as it follows

\begin{align*} \int_{-\infty}^0\frac{xe^x\operatorname{gd}(x)}{1+e^{2x}}\mathrm dx&=\int_0^{\infty}\frac{xe^x\operatorname{gd}(x)}{1+e^{2x}}\mathrm dx\\ &=\int_0^{\infty}\frac{xe^x\left[2\arctan(e^x)-\frac\pi2\right]}{1+e^{2x}}\mathrm dx\\ &=\int_0^{\infty}\frac{xe^{-x}\left[\frac\pi2-2\arctan(e^{-x})\right]}{1+e^{-2x}}\mathrm dx\\ &=-\int_1^0\frac{\log(y)y\left[\frac\pi2-2\arctan(y)\right]}{1+y^2}\frac{\mathrm dy}y\\ &=\int_0^1\frac{\log(y)\left[2\arctan(y)-\frac\pi2\right]}{1+y^2}\mathrm dy\\ \end{align*}

Now we are left with a rather "simple" integral. Enforcing the substitution $\arctan(y)\mapsto y$ this further reduces to

\begin{align*} \int_0^1\frac{\log(y)\left[2\arctan(y)-\frac\pi2\right]}{1+y^2}\mathrm dy&=\int_0^{\pi/4}\log(\tan y)\left[2y-\frac\pi2\right]\mathrm dy\\ &=2\int_0^{\pi/4}\left[y-\frac\pi4\right]\log(\tan y)\mathrm dy \end{align*}

Applying Integration By Parts with $u=y-\frac\pi4$ and $\mathrm dv=\log(\tan y)$. For this purpose we may introduce the Clausen Function $\operatorname{Cl}_2(z)$ which is capable of expressing the anti-derivative of $\log(\tan y)$. Thus, we get

\begin{align*} 2\int_0^{\pi/4}\left[y-\frac\pi4\right]\log(\tan y)\mathrm dy&=2\underbrace{\left[\left(y-\frac\pi4\right)\left(-\frac12\operatorname{Cl}_2(2y)-\frac12\operatorname{Cl}_2(\pi-2y)\right)\right]_0^{\pi/4}}_{=0}\\ &~+\int_0^{\pi/4}\operatorname{Cl}_2(2y)+\operatorname{Cl}_2(\pi-2y)\mathrm dy\\ &=\frac12\int_0^{\pi/2}\operatorname{Cl}_2(y)\mathrm dy+\frac12\int_0^{\pi/2}\operatorname{Cl}_2(\pi-y)\mathrm dy\\ &=\frac12\int_0^{\pi/2}\operatorname{Cl}_2(y)\mathrm dy+\frac12\int_{\pi/2}^\pi\operatorname{Cl}_2(y)\mathrm dy\\ &=\frac12\int_0^\pi\operatorname{Cl}_2(y)\mathrm dy\\ &=\frac12[\zeta(3)-\operatorname{Cl}_3(\pi)]\\ &=\frac12[\zeta(3)+\eta(3)]\\ &=\frac12\left[\zeta(3)+\frac34\zeta(3)\right]\\ &=\frac78\zeta(3) \end{align*}

$$\therefore~\int_{-\infty}^0\frac{xe^x\operatorname{gd}(x)}{1+e^{2x}}\mathrm dx~=~\frac78\zeta(3)$$

The Clausen Function is just a more convenient way to deal with the Fourier Series Expansion of $\log(\tan y)$ which would lead to the same result. The here used properties of the Clausen Function are quite easy to be proven by using the series representation and integrating termwise within the given borders of integration.