<s>[INST] Using complex analysis, show that $\int_0^\infty\frac{\arctan(x)}{1+x^2}\,dx=\frac{\pi^2}8$ The result is of course confirmed by substituting $x\mapsto\tan(x)$:
$$I = \int_0^\infty \frac{\arctan(x)}{1+x^2} \, dx = \int_0^{\frac\pi2} x \, dx = \frac{\pi^2}8$$
Or integrating by parts:
$$I = \lim_{x\to\infty} \arctan^2(x) - I \implies 2I = \left(\frac\pi2\right)^2 \implies I = \frac{\pi^2}8$$
Or splitting the integral at $x=1$ and substituting $x\mapsto\frac1x$ on the integral over $[1,\infty)$:
$$I= \int_0^1 \frac{\arctan(x) + \arctan\left(\frac1x\right)}{1 + x^2} \, dx = \frac\pi2 \int_0^1 \frac{dx}{1+x^2} = \frac{\pi^2}8$$
Or differentiating under the integral sign:
$$I(a) = \int_0^\infty \frac{\arctan(ax)}{1+x^2} \, dx \implies I'(a) = \int_0^\infty \frac x{(1+x^2)(1+a^2x^2)} \, dx = \frac{\ln(a)}{a^2-1} \\ I(0) = 0 \implies I(1) = \int_0^1 \frac{\ln(x)}{x^2-1} \, dx = \frac{\pi^2}8$$
Or getting the same integral of $\frac{\ln(x)}{x^2-1}$ by converting $\arctan(x)$ to an integral representation and computing the resulting double integral (per @Dr.WolfgangHintze's suggestion) using the same substitution as in the third method above:
$$\begin{align*}
I &= \int_0^\infty \int_0^x \frac x{(1+x^2)(1+x^2y^2)} \, dy \, dx \\[1ex]
&= \int_0^\infty \int_y^\infty \frac x{(1+x^2)(1+x^2y^2)} \, dx \, dy \\[1ex]
&= \frac12 \int_0^\infty \frac{\ln\left(\frac{y^2+y^4}{1+y^4}\right)}{y^2-1} \, dy \\[1ex]
&= \int_0^\infty \frac{\ln(y)}{y^2-1} \, dy + \frac12 \int_0^\infty \frac{\ln(1 + y^2)}{y^2-1} \, dy - \frac12 \int_0^\infty \frac{\ln(1+y^4)}{y^2-1} \, dy \\[1ex]
&= (1 + 2 - 2) \int_0^1 \frac{\ln(y)}{y^2-1} \, dy = \frac{\pi^2}8
\end{align*}$$

I was wondering how, if at all possible, one might approach it with the residue theorem? I see that $z=\pm i$ are simple poles of $\frac1{1+z^2}$, but they're also the branch points of $\arctan(z)$, since
$$\arctan(z) = -\frac i2 \log\left(\frac{i-z}{i+z}\right)$$
so I don't believe the theorem can be readily applied here. The integrand is odd so I don't think there's much to infer from symmetry. Maybe there's a way to massage the integrand to get closer to something with which we can use a contour integral.
 [/INST] Complex integration can also be used to evaluate the integral; often, this is a very power tool for finding general solutions - if the system has an appropriate symmetry.
We will consider a general case:
$$I(a,b)=\int_0^\infty\frac{x^b\arctan x}{a^2+x^2}dx;\,b\in[0;1);\,a>0\tag{0}$$
We will also choose $a<1$ - it allows us to separate poles from branch points; as the integrand is a continuous function of $a$, our solution will be valid for all $a>0$ (after an appropriate rearrangement if $a>1$).
We define $\arctan z$ in a standart way: $\displaystyle\arctan z=-\frac{i}{2}\ln\frac{i-z}{i+z}$, and the regular cuts are $[-i;-i\infty);\,[i;i\infty)$.
The modified keyhole contour looks

To define the integrand on the banks of the cut (upper half-plane), we choose $z=r<1$ on the axis $X$, then make a turn counter-clockwise around $\displaystyle z=0$: $\displaystyle\,r\to re ^\frac{\pi i}{2}; \ln z\to-\frac{i}{2}\ln\frac{i-ir}{i+ir}=-\frac{i}{2}\ln\frac{1-r}{1+r}$. To move on the bank of the cut (area A), we make another turn (angle $\pi$ around $z=i$; counter-clockwise): $\displaystyle 1-r\to (1-r)e ^{\pi i};\,-\frac{i}{2}\ln\frac{1-r}{1+r}\to -\frac{i}{2}\ln\Big(\frac{1-r}{1+r}e ^{\pi i}\Big)=-\frac{i}{2}\ln\frac{1-r}{1+r}+\frac{\pi}{2}$
If we move to the opposite bank (region C), we will get $-\frac{i}{2}\ln\frac{1-r}{1+r}-\frac{\pi}{2}$. On the both banks of the cut $z= re ^\frac{\pi i}{2}$ and $z^b= re^\frac{\pi ib}{2}$.
Integral along our close contour
$$\oint=I(a,b)(1-e^{2\pi ib})+I_A+I_C+I_D+I_E+I_R+I_{r_k}$$
It is not difficult to show that the integrals along a big circle and small circles (around $z=0,\pm i$) $I_R$ and $I_{r_k}$ tend to zero at $R\to\infty$ and $r_k\to0$. Making the similar analysis for the bottom part of the contour (as was done above for the regions A, C) and bearing in mind that in this case $z= re ^\frac{3\pi i}{2}$ (our contour goes counter-clockwise)
$$I(a,b)(1-e^{2\pi ib})+I_A+I_C+I_D+I_E=2\pi i\underset {z=ae ^\frac{\pi i}{2}, \,ae ^\frac{3\pi i}{2}}{\operatorname{Res}}-\frac{i}{2}\frac{z^b}{z^2+a^2}\ln\frac{i-z}{i+z}\tag{1}$$
To evaluate $\displaystyle I_{A, C, D, T}$ we notice that the terms with $\ln\frac{i-z}{i+z}$ cancel (we integrate in the opposite directions), and we are left, for example
$$I_A+I_C=-\frac{\pi}{2}e ^\frac{\pi i}{2}e ^\frac{\pi ib}{2}\int_1^\infty\frac{x^b}{a^2-x^2}dx+\frac{\pi}{2}e ^\frac{\pi i}{2}e ^\frac{\pi ib}{2}\int_\infty^1\frac{x^b}{a^2-x^2}dx$$
$$=-\pi ie ^\frac{\pi ib}{2}\int_1^\infty\frac{x^b}{a^2-x^2}dx$$
In the similar way,
$$I_D+I_E=-\pi ie ^\frac{3\pi ib}{2}\int_1^\infty\frac{x^b}{a^2-x^2}dx$$
The residues evaluation gives
$$2\pi i\sum\operatorname {Res}=-\frac{\pi i}{2}a^{b-1}\ln\frac{1-a}{1+a}\Big(e ^\frac{\pi ib}{2}+e ^\frac{3\pi ib}{2}\Big)$$
Putting all in (1)
$$I(a,b)=\frac{\pi}{2\sin\frac{\pi b}{2}}\int_1^\infty\frac{x^b}{x^2-a^2}dx-\frac{\pi}{4\sin\frac{\pi b}{2}}a^{b-1}\ln\frac{1+a}{1-a}$$
Integrating the first term by part
$$\boxed{\,\,I(a,b)=\frac{\pi}{4a\sin\frac{\pi b}{2}}\bigg((1-a^b)\ln\frac{1+a}{1-a}+b\int_1^\infty x^{b-1}\ln\frac{x+a}{x-a}dx\bigg);\, a\in(0;1]\,\,}\tag{2}$$
Now we can consider different cases.

*

*Leading $\displaystyle b\to 0$
$$I(a,0)=\int_0^\infty\frac{\arctan x}{a^2+x^2}dx=\int_1^\infty\frac{\ln x}{x^2-a^2}dx-\frac{\ln a}{2a}\ln\frac{1+a}{1-a}\tag{a}$$

*Putting in turn $\displaystyle a=1$
$$I(1,0)=\int_0^\infty\frac{\arctan x}{1+x^2}dx=\int_1^\infty\frac{\ln x}{x^2-1}dx=\frac{\pi^2}{8}\tag{b}$$

*Taking $a\to1$ (the first term in (2) is zero), integrating the second term by part, using expressions for gamma (duplication formula) and digamma functions, we can get an interesting result
$$I(1,b)=\int_0^\infty\frac{x^b\arctan x}{1+x^2}dx=\frac{\pi}{4\sin\frac{\pi b}{2}}\bigg(\psi\Big(\frac{1}{2}\Big)-\psi\Big(\frac{1-b}{2}\Big)\bigg);\,b\in[0;1)\tag{c}$$
When $b\to 0$, we get again
$$I(1;0)=\int_0^\infty\frac{\arctan x}{1+x^2}=\frac{1}{4}\psi'\Big(\frac{1}{2}\Big)=\sum_{k=1}^\infty\frac{1}{(2k-1)^2}=\frac{\pi^2}{8}\tag{d}$$

*The integral has also a closed form for $\displaystyle b=\frac{1}{2}$ (and arbitrary $a>0$):
$$I\Big(a, \frac{1}{2}\Big)=\frac{\pi}{2\sqrt2 a}\bigg((1-\sqrt a)\ln\frac{1+a}{1-a}+\frac{1}{2}\int_1^\infty\ln\frac{x+a}{x-a}\frac{dx}{\sqrt x}\bigg)$$
Making in the second term the substitution $t=\sqrt x$, integrating by part and using $\arctan\frac{1}{s}=\frac{\pi}{2}-\arctan s$
$$= \frac{\pi}{2\sqrt2 a}\bigg((1-\sqrt a)\ln\frac{1+a}{1-a}+\ln\frac{1-a}{1+a}+\sqrt a\ln\frac{1+\sqrt a}{1-\sqrt a}+2\sqrt a\arctan\sqrt a\bigg)$$
and we get the known result:
$$I\Big(a, \frac{1}{2}\Big)=\int_0^\infty\frac{\sqrt x\arctan x}{a^2+x^2}dx=\frac{\pi}{\sqrt {2a}}\bigg(\arctan\sqrt a+\ln\frac{1+\sqrt a}{\sqrt{1+a}}\,\bigg),\,a>0\tag{e}$$
</s>