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.
Best Answer
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.