Calculating $\int_0^{\infty } \left(\text{Li}_2\left(-\frac{1}{x^2}\right)\right)^2 \, dx$

With a substitution and a step of integration by parts, the problems boils down to computing: $$ I = -4\int_{0}^{+\infty}\frac{\log(1+x^2)\,\text{Li}_2(-x^2)}{x^2}\,dx. \tag{1}$$ Integrating by parts again, the problem boils down to computing: $$ I_1 = \int_{0}^{+\infty}\frac{\log(1+x^2)^2}{x^2}\,dx,\qquad I_2=\int_{0}^{+\infty}\arctan\left(\frac{1}{x}\right)\frac{\log(1+x^2)}{x}\,dx.\tag{2} $$ The first integral is straightforward: $$ I_1 = \frac{1}{2}\left.\frac{d^2}{d\alpha^2}\int_{0}^{+\infty}\frac{(1+x)^{\alpha}-1}{x^{3/2}}\,dx\,\right|_{\alpha=0} \tag{3}$$ since the innermost integral can be evaluated in terms of the beta function.

That leads to $I_1=4\pi\log 2$. Now we just need to compute: $$ J = \int_{1}^{+\infty}\int_{0}^{+\infty}\frac{\log(1+x^2)}{1+t^2 x^2}\,dx\,dt \tag{4}$$ to recover the value of $I_2=\frac{\pi^3}{8}$. That proves the conjecture:

$$ I = \frac{4\pi^3}{3}+32\pi\log 2.\tag{5}$$

Another way to evaluate the integral $$I_{2} = \int_{0}^{\infty} \frac{\arctan (\frac{1}{x}) \log(1+x^{2})}{x} \, dx$$ in Jack D'Aurizio's answer (which according to Wolfram Alpha evaluates to $\frac{\pi^{3}}{12}$) is to consider the complex function $$f(z) = \frac{\arctan \left(\frac{1}{z} \right) \log(1-iz)}{z} = \frac{\text{arccot}(z) \log(1-iz)}{z}.$$

Using the principal branch of the logarithm, $\text{arccot}(z)$ has a branch cut on $[-i, i]$ and $\log(1-iz)$ has a branch cut on $(-i\infty, -i]$.

So integrating around a semicircle in the upper half-plane deformed around the branch cut and using the fact $\text{arccot}(z) \sim \frac{1}{z}$ when $z$ is large in magnitude, we get

$$\int_{-\infty}^{\infty} \frac{\arctan(\frac{1}{x}) \log(1-ix)}{x} \, dx + \int_{0}^{1} \frac{\frac{i}{2} (2 \pi i) \log(1+t)}{it} \, i \, dt =0.$$

Then equating the real parts on both sides of the equation, we get

$$\int_{0}^{\infty} \frac{\arctan(\frac{1}{x}) \log(1+x^{2})}{x} \, dx = \pi \int_{0}^{1} \frac{\log(1+t)}{t} \, dt = - \pi \, \text{Li}_{2}(-1) = \frac{\pi^{3}}{12}.$$

Ok i want to tackle the problem from a different starting point:

Recall that the $\text{Li}_2(x)$ has the following integral representation which is closely related to the Debye functions, well known in condensed matter physics:

$$ \text{Li}_2(z)=\frac{z}{\Gamma(2)}\int_{0}^{\infty}\frac{t}{e^t-z}dt \quad \quad(1) $$

This leads us to the following representation of our Integral $$ I=\int_{0}^{\infty}\left(\text{Li}_2\left(\frac{1}{x^2}\right)\right)^2dx=\\ \int_{0}^{\infty}\mathrm{d}t_1t_1\int_{0}^{\infty}\mathrm{d}t_2t_2\int_{0}^{\infty}dx\frac{1}{x^2 e^{t_1}+1}\frac{1}{x^2 e^{t_2}+1} $$

The inner integral is quite a straightforward exercise in residue calculus and yields (using parity)

$$ I=8\pi\int_{0}^{\infty}\mathrm{d}t_1t_1\int_{0}^{\infty}\mathrm{d}t_2t_2\frac{1}{e^{t_1}+e^{t_2}} $$

Using (1) again to perform for example the $t_1$ integral, we get

$$ I=8\pi\int_{0}^{\infty}\mathrm{d}t_2t_2 e^{-t_2} \text{Li}_2\left(-e^{t_2}\right) $$

Defining $p(t)=t e^{-t} $ and $ q(t)=\text{Li}_2\left(-e^{t}\right)$ (relabeling $t=t_2$) we perform two integration by parts ending up with ( $p^{(-n)}(x)$ denotes the $n$ primitive with respect to $x$)

$$ \frac{1}{8\pi}I=\underbrace{[p^{(-1)}(t)q(t)-p^{(-2)}(t)q'(t)]_0^{\infty}}_{=\frac{\pi ^2}{12}+2\log (2)}+\underbrace{\int_0^{\infty}\overbrace{\frac{t+2}{e^t+1}}^{p^{(-2)}(t)q"(t)}dt}_{J} $$

Here we have used the special value $\text{Li}_2\left(-1\right)=-\frac{\pi^2}{12}\quad (2)$ (See also the appendix for a proof of this fact). The $J$ is also straightforwardly calculated in the appendix and yields

$$ J=\frac{\pi ^2}{12}+2\log (2) $$

Ande therefore

$$ I=\frac{4\pi^3}{3}+32\pi \log(2) $$

Q.E.D

Appendix

Calculation of $J$

$$ J=2\underbrace{\int_{0}^{\infty}\frac{1}{1+e^t}}_{J_1}+\underbrace{\int_{0}^{\infty}\frac{t}{1+e^t}}_{J_2} $$

$J_1$ can be done by done by letting $e^t=z$

$$ J_1=\int_{1}^{\infty}\frac{1}{z(z+1)}=[\log (z)-\log (z+1)]_0^{\infty}=\log(2) $$

$J_2$ may be done by different methods, but most straightforwardly i think by using geometric series and integrating term wise

$$ J_2 =\sum_{n=1}^{\infty}(-1)^n\int_{0}^{\infty}e^{-(n+1)x}=\sum_{n=1}^{\infty}\frac{(-1)^{n+1}}{n^2}=\eta(2)=\frac{\pi^2}{12} \quad (3) $$

here we used a special value of the Dirichlet $\eta$-function

Note that by a change of variables $e^{-z}=q$ the last integral in (3) we also have $\sum_{n=1}^{\infty}\frac{(-1)^{(n+1)}}{n^2}=-\text{Li}_2(-1)$ which proves (2).

Using $J=2J_1+J_2$ the stated result follows