Can $\int_0^{\infty} \frac{\ln ^n x}{\left(1-x^2\right)^2} d x$ be expressed in terms of the zeta function of an even integer? 
Being attracted by the exact values of the integrals, $$
\begin{aligned}
& \int_0^{\infty} \frac{\ln ^2 x}{\left(1-x^2\right)^2} d x=\frac{\pi^2}{4} \\
& \int_0^{\infty} \frac{\ln ^3 x}{\left(1-x^2\right)^2} d x=-\frac{\pi^4}{16} \\
& \int_0^{\infty} \frac{\ln ^4 x}{\left(1-x^2\right)^2}dx=\frac{\pi^4}{4} \\
& \int_0^{\infty} \frac{\ln ^5 x}{\left(1-x^2\right)^2}dx=-\frac{\pi^6}{8} \\
& \int_0^{\infty} \frac{\ln ^6 x}{\left(1-x^2\right)^2}dx=\frac{3 \pi^6}{4},
\end{aligned}
$$
I dare guess that the integral in general
$$I_n=\int_0^{\infty} \frac{\ln ^n x}{\left(1-x^2\right)^2}dx$$
may be expressed in terms of the zeta function of an even integers.


As usual, we split it into two integrals as
$$
\int_0^{\infty} \frac{\ln ^n x}{\left(1-x^2\right)^2} d x=\int_0^1 \frac{\ln ^n x}{\left(1-x^2\right)^2} d x+\int_1^{\infty} \frac{\ln ^n x}{\left(1-x^2\right)^2} d x
$$
$$\because \int_1^{\infty} \frac{\ln ^n x}{\left(1-x^2\right)^2} d x \stackrel{x\mapsto\frac{1}{x}}{=} (-1)^n \int_0^1 \frac{x^2 \ln ^nx}{\left(1-x^2\right)^2} d x \\= (-1)^n \int_0^1 \frac{\ln ^n x}{\left(1-x^2\right)^2} d x+(-1)^{n+1} \int_0^1 \frac{\ln ^n x}{1-x^2} d x $$
$$ \therefore \int_0^{\infty} \frac{\ln ^n x}{\left(1-x^2\right)^2} d x=(1+(-1)^n) \int_0^1 \frac{\ln ^n x}{\left(1-x^2\right)^2} d x +(-1)^{n+1} \int_0^1 \frac{\ln ^n x}{1-x^2} d x  $$

So when $n$ is odd, $$
\int_0^{\infty} \frac{\ln ^n x}{\left(1-x^2\right)^2} d x=\int_0^1 \frac{\ln ^n x}{1-x^2} dx= -n !\left(1-\frac{1}{2^{n+1}}\right) \zeta(n+1)
$$
where the last answer comes from my post 1.

When $n$ is even, $$
\begin{aligned}
\int_0^{\infty} \frac{\ln ^n x}{\left(1-x^2\right)^2}d x&=2 \int_0^1 \frac{\ln ^n x}{\left(1-x^2\right)^2} d x-\int_0^1 \frac{\ln ^n x}{1-x^2} d x \\
& =n !\left[\left(1-\frac{1}{2^n}\right) \zeta(n)+\left(1-\frac{1}{2^{n+1}}\right) \zeta(n+1)\right]-n !\left(1-\frac{1}{2^{n+1}}\right)\zeta(n+1) \\
& =n !\left(1-\frac{1}{2^n}\right) \zeta(n)
\end{aligned}
$$
where the second last line comes from my post 2 and post 1.

We can now conclude that
$$
\int_0^{\infty} \frac{\ln ^n x}{\left(1-x^2\right)^2} d x=\left\{\begin{array}{l}
-n !\left(1-\frac{1}{2^{n+1}}\right) \zeta(n+1) \quad  \text { if }n \textrm{  is odd greater than }1. \\
n !\left(1-\frac{1}{2^n}\right) \zeta(n) \qquad\qquad    \text { if } n \textrm{  is even. }
\end{array}\right.
$$
can be expressed in terms of the zeta function of even integer.
Your comments and alternative solutions are highly appreciated.
 A: $\newcommand{\d}{\,\mathrm{d}}\newcommand{\J}{\mathcal{J}}\newcommand{\H}{\mathcal{H}}\newcommand{\pv}{\operatorname{P.V.}}$Here is one alternative solution, using complex analysis.
Let's remove the nasty double pole from the denominator: $$\begin{align}I_n&=\int_0^\infty\frac{\ln^nx}{(1-x^2)^2}\\&=(-1)^n\int_0^\infty\frac{\ln^nx}{(1-x^2)^2}\cdot x^2\d x\\&=(-1)^n2^{-(n+1)}\int_0^\infty\frac{\ln^nx}{(x-1)^2}\sqrt{x}\d x\\&=(-1)^n2^{-(n+1)}\int_0^\infty\frac{1}{x-1}\left(\frac{n\ln^{n-1}x}{x}\sqrt{x}+\frac{1}{2}\frac{\ln^nx}{\sqrt{x}}\right)\d x\\&=(-1)^n2^{-(n+2)}(2n\cdot M_{n-1}+M_n)\end{align}$$Where integration by parts on $(x-1)^{-2}\to-(x-1)^{-1}$ was used and: $$M_k:=\int_0^\infty\frac{\ln^kx}{(x-1)\sqrt{x}}\d x$$For integer $k\ge1$.
Fix $s\in\Bbb C$ with $-\frac{1}{2}<\Re s<\frac{1}{2}$. We consider:
$$\J(s):=\int_0^\infty\frac{x^s\log x}{(x-1)\sqrt{x}}\d x$$
Where both $\log$ and $x^s$ are computed via, $0\le\arg z<2\pi$. Let $f(z)=\frac{z^{s-1/2}\log z}{z-1}$ under these definitions.
Two important facts: $\lim_{|z|\to\infty}|z\cdot f(z)|=0$ and $\lim_{|z|\to0^+}|z\cdot f(z)|=0$. Consider:
$$\lim_{R\to\infty}\oint_{\H_R}f(z)\d z$$
Where $\H_R$ is the truncated Hankel contour: for some fixed $\epsilon>0$, we take $\H_R$ to run $R-i\epsilon\to-i\epsilon$ in a straight line, $-i\epsilon\to i\epsilon$ in a semicircular arc that wraps around zero (it passes through $-\epsilon$, e.g.) and goes $i\epsilon\to R+i\epsilon$ in a straight line again. Since $f$ is holomorphic on $\Bbb C\setminus[0,\infty)$. We can deform $\H_R$ to a very large circle of radius $\sqrt{R^2+\epsilon^2}$ without crossing $[0,\infty)$ (keeping the endpoints $R\pm i\epsilon$ fixed). As $f$ is holomorphic away from that line, this does not change the value of the integral. It is then clear that $\oint_{H_R}f\to0$.
But let's compute this limit another way. Fix $\delta\in(0,1)$.  As $R$ grows larger we can push the top horizontal line down to the real axis, squeeze the semicircle down to a point, and then push the bottom line down so it looks like $\infty\to1+\delta$ in a straight line, $1+\delta\to1-\delta$ in a semicircular arc dipping below the real axis, and then $1-\delta\to0$ in a straight line. Now, the integral over the semicircle at zero vanishes as we squish the semicircle. The integral on the top line, as it descends to the real axis, is precisely $\J(s)$. Being careful with branch jumps, we see the integral on the bottom line tends to:
$$-\left(\int_0^{1-\delta}+\int_{1+\delta}^\infty\right)\frac{x^{s-1/2}\cdot e^{2\pi i(s-1/2)}\cdot(\log x+2\pi i)}{x-1}\d x-\oint_{C_\delta}\frac{z^{s-1/2}\log z}{z-1}\d z$$
The important point is that these (limiting) deformations of the contour took place within the holomorphic domain of $f$, so did not change the value of the limit $\lim_R\oint_{\H_R}f$. Moreover, we get an equality for any $\delta\in(0,1)$, so it is valid to take $\delta\to0^+$ without changing the equality.
Using the fact that $z^{s-1/2}\log z=-2\pi i\cdot e^{2\pi is}+\mathcal{O}(z-1)$ as $z\to1$ from below) we see, as $\delta\to0^+$, the contribution from the lower line of the contour can be taken to be:
$$e^{2\pi is}\J(s)+2\pi i\cdot e^{2\pi is}\pv\int_0^\infty\frac{x^{s-1/2}}{x-1}\d x-2\pi^2e^{2\pi is}$$Using the fact that the sum of all contributions is zero, we deduce: $$\begin{align}(1+e^{2\pi is})\J(s)&=2\pi\cdot e^{2\pi is}(\pi-i\pv)\\\J(s)&=\pi\cdot e^{\pi is}\sec(\pi s)\cdot(\pi-i\pv)\\&=\pi(1+i\tan\pi s)(\pi-i\pv)\end{align}$$
Since $\J(s)$ must be real when $s$ is real, we see it is actually forced that the principal value is $\pi\tan\pi s$ for real $s$. Then by analytic continuation, for all $-\frac{1}{2}<\Re s<\frac{1}{2}$ we get the value of $\pi\tan\pi s$. This principal value can also be calculated with complex analysis, but it’s nice to get it ‘for free’ like this.
This concludes: $$\J(s)=\pi^2\sec^2\pi s=\pi^2+2\pi^2\sum_{n\ge1}(-1)^n2^{2n}(2^{2(n+1)}-1)B_{2(n+1)}\frac{(\pi s)^{2n}}{(n+1)(2n)!}$$
Differentiating and evaluating at $s=0$ allows us to read off the values of $M_k$ as: $M_{2k}=0$ for all $k$, and: $$M_{2n-1}=\zeta(2n)\frac{(2n)!(2^{2n}-1)}{n}$$
For all integer $n\ge0$, using Euler's relation between the Bernoulli numbers and $\zeta$.
Let $n=2m-1$. Using the formulae at the top of the post, we find: $$I_{2m-1}=-2^{-(2m+1)}M_{2m-1}=-(2m-1)!(1-2^{-2m})\cdot\zeta(2m)$$
If $n=2m$, we similarly find: $$I_{2m}=(2m!)(1-2^{-2m})\cdot\zeta(2m)$$
These agree with your formula.
