Conjecture regarding integrals of the form $\int_0^\infty \frac{(\log{x})^n}{1+x^2}\,\mathrm{d}x$. I have been playing around a bit with integrals of the form $$I(n)=\int_0^\infty \frac{(\log{x})^n}{1+x^2}\,\mathrm{d}x,\,\,n\in\mathbb{Z}^+,$$ and I am trying to obtain a closed form solution for $I(n).$ I believe the special cases $I(1)$ and $I(2)$ are somewhat well-known, but I will go over them. When $n=1,$ we have $$I(1)=\int_0^\infty \frac{\log{x}}{1+x^2}\,\mathrm{d}x=\int_0^1 \frac{\log{x}}{1+x^2}\,\mathrm{d}x+\int_1^\infty \frac{\log{x}}{1+x^2}\,\mathrm{d}x.$$ This can be easily shown to be zero by performing the substitution $x=1/y,$ which will yield $$\int_0^1 \frac{\log{x}}{1+x^2}\,\mathrm{d}x=-\int_1^\infty \frac{\log{x}}{1+x^2}\,\mathrm{d}x.$$ Thus, $I(1)=0.$ Clearly, this can be generalized to all odd integers, and $I(2n+1)=0.$ In the case of $n=2$, first observe through the same substitution as above that $$\int_0^1 \frac{(\log{x})^2}{1+x^2}\,\mathrm{d}x=\int_1^\infty \frac{(\log{x})^2}{1+x^2}\,\mathrm{d}x.$$ This implies that $$I(2)=2\int_0^1 \frac{(\log{x})^2}{1+x^2}\,\mathrm{d}x,$$ which is an easier integral to work with. Performing the substitution $x=e^{-y}$ yields $$I(2)=2\int_0^\infty y^2\left(e^{-y}-e^{-3y}+e^{-5y}-\cdot\cdot\cdot\right)\,\mathrm{d}y.$$ Using the identity $$\int_0^\infty x^2 e^{-ax}=\frac{2}{a^3},$$ we obtain $$I(2)=4\left(\frac{1}{1^3}-\frac{1}{3^3}+\frac{1}{5^3}-\cdot\cdot\cdot\right)=4\cdot\frac{\pi^3}{32}=\frac{\pi^3}{8}.$$
I'm not sure how this infinite series is evaluated, but I found this result in a book. I used Mathematica to check a few more values, and I found that $I(4)=5\pi/32,\,I(6)=61\pi/128,\,$ and $I(8)=1385\pi/512.$ Clearly the pattern is $$I(2n)=A_{2n}\left(\frac{\pi}{2}\right)^{2n+1},$$ where $A_{2n}$ is some constant. It turns out that these constants are the Euler numbers, which are the coefficients $E_k$ corresponding to the series $$\operatorname{sech}x=\sum_{k=0}^\infty\frac{E_k}{k!}x^k.$$
All Euler numbers corresponding to odd $n$ are zero, and the first few even Euler numbers are $E_0=1,\, E_2=-1,\, E_4=5, \,E_6=-61,\,$ and $E_8=1385.$ Thus, I have conjectured that $$I(2n)=(-1)^n E_{2n} \left(\frac{\pi}{2}\right)^{2n+1}$$ for all $n\in\mathbb{Z}^+.$ I suppose this could be extended to $n\in\mathbb{Z}_{\geq 0}$ thusly:
$$I(n)=i^n E_n \left(\frac{\pi}{2}\right)^{n+1},$$ since $E_n=0$ for odd $n.$ So the question, of course, is how to prove this. I tried generalizing the method I used for $I(2),$ and I found that $$\int_0^\infty x^n e^{-ax}=\frac{n!}{a^{n+1}},\,\,n\in\mathbb{Z}_{\geq 0}.$$ Using this, I obtained $$I(2n)=n!\left(\frac{1}{1^{n+1}}-\frac{1}{3^{n+1}}+\frac{1}{5^{n+1}}-\cdot\cdot\cdot\right)=n!\sum_{k=0}^\infty\frac{(-1)^k}{(2k+1)^{n+1}}.$$ Mathematica wasn't able to evaluate this sum, even for the case of $n=2.$ It gives some expression involving multiple zeta functions, with which I have no experience. Even if we can't prove this, I would be interested to know why the Euler numbers might appear here. Any help would be greatly appreciated.
Edit:
As Claude Leibovici helped point out, there final series expression should be $$I(2n)=2(2n)!\sum_{k=0}^\infty \frac{(-1)^k}{(2k+1)^{2n+1}}=\frac{(2n)!}{2^{4n+1}}\left[\zeta\left(2n+1, \frac{1}{4}\right)-\zeta\left(2n+1, \frac{3}{4}\right)\right].$$

Integrate $f(z)=\dfrac{z^{s}}{1+z^2}$, with the branch cut placed on the positive real axis and $-1<\operatorname{Re}(s)<1$, along a keyhole contour deformed around $[0,R]$. Along the big arc of radius $R$,
$$0\leq\left|\ \int_{R\exp\left(i[0,2\pi]\right)} f(z)\ dz\ \right|\leq\frac{2\pi R^{\operatorname{Re}(s)+1}}{R^2-1}\to0 \text{ as }R\to\infty$$
and along the small arc of radius $\epsilon$,
$$0\leq\left|\ \int_{\epsilon\exp\left(i[0,2\pi]\right)} f(z)\ dz\ \right|\leq\frac{2\pi \epsilon^{\operatorname{Re}(s)+1}}{1-\epsilon^2}\to0 \text{ as }\epsilon\to 0$$
so taking $R\to\infty$ and $\epsilon\to 0$ and applying the Residue Theorem,
\begin{align}
(1-e^{2\pi i s})\int^\infty_0\frac{x^{s}}{1+x^2}\ dx
&=\pi\left(e^{\pi is/2}-e^{3\pi i s/2}\right)\\
\implies \int^\infty_0\frac{x^{s}}{1+x^2}\ dx
&=\pi\cdot\frac{e^{\pi is/2}-e^{-\pi i s/2}}{e^{\pi is}-e^{-\pi is}}
=\pi\cdot\frac{\sin\left(\frac{\pi s}{2}\right)}{\sin\left(\pi s\right)}
=\frac{\pi}{2}\sec\left(\frac{\pi s}{2}\right)
\end{align}
Therefore,
\begin{align}
I(2n)
&=\left.\frac{\pi}{2}\frac{d^{2n}}{ds^{2n}}\sec\left(\frac{\pi s}{2}\right)\right|_{s\to 0}\\
&=\left.\frac{\pi}{2}\frac{d^{2n}}{ds^{2n}}\sum^\infty_{k=0}\frac{(-1)^kE_{2k}}{(2k)!}\left(\frac{\pi s}{2}\right)^{2k}\right|_{s\to 0}\\
&=\left.\frac{\pi}{2}(2n)![s^{2n}]\sum^\infty_{k=0}\frac{(-1)^kE_{2k}}{(2k)!}\left(\frac{\pi s}{2}\right)^{2k}\right|_{s\to 0}\\
&=\frac{\pi}{2}(2n)!\frac{(-1)^nE_{2n}}{(2n)!}\left(\frac{\pi}{2}\right)^{2n}\\
&=(-1)^nE_{2n}\left(\frac{\pi}{2}\right)^{2n+1}
\end{align}
