Ok quiero abordar el problema desde un punto de partida diferente:
Recordemos que el $\text{Li}_2(x)$ tiene la siguiente representación integral
que está estrechamente relacionado con el de Debye funciones, bien conocido en física de la materia condensada:
$$
\text{Li}_2(z)=\frac{z}{\Gamma(2)}\int_{0}^{\infty}\frac{t}{e^t-z}dt \quad \quad(1)
$$
Esto nos lleva a la siguiente representación de nuestros 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}
$$
El interior de la integral es bastante sencillo ejercicio en el residuo de cálculo
y de los rendimientos (utilizando la paridad)
$$
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}}
$$
El uso de (1) de nuevo para realizar, por ejemplo, el $t_1$ integral, obtenemos
$$
I=8\pi\int_{0}^{\infty}\mathrm{d}t_2t_2 e^{-t_2} \text{Li}_2\left(-e^{t_2}\right)
$$
La definición de $p(t)=t e^{-t} $ $ q(t)=\text{Li}_2\left(-e^{t}\right)$ (reetiquetar $t=t_2$)
se realizan dos de la integración por partes para terminar con ( $p^{(-n)}(x)$ indica el $n$ primitiva con respecto a $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}
$$
Aquí hemos utilizado el especial valor de $\text{Li}_2\left(-1\right)=-\frac{\pi^2}{12}\quad (2)$ (Véase también el apéndice para una prueba de este hecho).
El $J$ es también directamente calculado en el apéndice y los rendimientos
$$
J=\frac{\pi ^2}{12}+2\log (2)
$$
Ande, por tanto,
$$
I=\frac{4\pi^3}{3}+32\pi \log(2)
$$
Q. E. D
Apéndice
Cálculo de $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$ se logra dejando $e^t=z$
$$
J_1=\int_{1}^{\infty}\frac{1}{z(z+1)}=[\log (z)-\log (z+1)]_0^{\infty}=\log(2)
$$
$J_2$ puede realizar por diferentes métodos, pero la mayoría de principio a fin, creo que por el uso de series geométricas y la integración de plazo sabio
$$
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)
$$
aquí hemos utilizado un valor especial de la Dirichlet $\eta$-función
Tenga en cuenta que por un cambio de las variables de $e^{-z}=q$ la última integral en (3) también contamos $\sum_{n=1}^{\infty}\frac{(-1)^{(n+1)}}{n^2}=-\text{Li}_2(-1)$, lo que demuestra (2).
El uso de $J=2J_1+J_2$ declaró el resultado de la siguiente manera