Esta no es una respuesta completa, sino sólo una descripción de dos ideas que podrían ayudar a la evaluación de la integral $$ I \equiv 4 \int \limits_0^1 \left(\frac{\operatorname{li}(x)}{x}\right)^3 (x-1) \, \mathrm{d} x \, . $$ Se basan en métodos que se pueden aplicar para encontrar la integral más fácil $$ J \equiv \int \limits_0^1 \left(\frac{\operatorname{li}(x)}{x}\right)^2 \, \mathrm{d} x \, . $$
El primer enfoque se basa en la integración por partes y la serie $$ x-1 = \sum \limits_{k=1}^\infty \frac{1}{k!} \ln^k (x) \, , \, x > 0 \, .$$ Para evaluar $J$ podemos utilizar la antiderivada $x \mapsto 1-\frac{1}{x}$ de $x \mapsto \frac{1}{x^2}$ para evitar problemas con la singularidad de $\operatorname{li}(x)$ en $x = 1$ . Obtenemos \begin{align} J &= 2 \int \limits_0^1 \frac{\operatorname{li}(x)}{x} \frac{1-x}{\ln(x)} \, \mathrm{d} x = - 2 \sum \limits_{k=1}^\infty \frac{1}{k!} \int \limits_0^1 \frac{\operatorname{li}(x)}{x} \ln^{k-1} (x) \, \mathrm{d} x\\ &= 2 \sum \limits_{k=1}^\infty \frac{1}{k! k} \int \limits_0^1 \ln^{k-1} (x) \, \mathrm{d} x = 2 \sum \limits_{k=1}^\infty \frac{1}{k! k} (-1)^{k-1} (k-1)! \\ &= 2 \sum \limits_{k=1}^\infty \frac{(-1)^{k-1}}{k^2} = 2 \eta (2) = \zeta(2) = \frac{\pi^2}{6} \, . \end{align} Del mismo modo, podemos utilizar la antiderivada $x \mapsto \frac{(x-1)^2}{2 x^2}$ de $x \mapsto \frac{x-1}{x^3}$ para encontrar \begin{align} I &= - \frac{3}{2} \int \limits_0^1 \left(\frac{\operatorname{li}(x)}{x}\right)^2 \frac{(x-1)^2}{\ln(x)} \, \mathrm{d} x \\ &= \frac{3}{2} \sum \limits_{k=0}^\infty \frac{1}{k!} \int \limits_0^1 \operatorname{li}^2 (x) \frac{1-x}{x} \frac{\ln^{k-1} (x)}{x} \, \mathrm{d} x \, . \tag{A} \end{align} Ahora podemos integrar por partes una vez más para obtener al menos un término que reduzca a un múltiplo de $\zeta(3)$ como en el caso más sencillo. Sin embargo, aún no he conseguido resolver las integrales restantes. Por supuesto, podríamos volver a utilizar la serie para expresar las restantes $1-x$ en términos de potencias de logaritmos, pero eso no parece resolver el problema.
La segunda sugerencia emplea la serie de Fourier-Laguerre $$ \operatorname{li} (x) = - x \sum_{n=0}^\infty \frac{\mathrm{L}_n (-\ln(x))}{n+1} \, , \, x \in (0,1) \, , \tag{B}$$ de la integral logarítmica. Se puede demostrar derivando una relación de recurrencia para los coeficientes a partir de la de los polinomios de Laguerre.
Utilizando la sustitución $x = \mathrm{e}^{-t}$ y la relación de ortogonalidad de los polinomios de Laguerre obtenemos inmediatamente $$ J = \sum \limits_{p=0}^\infty \sum \limits_{q=0}^\infty \frac{1}{(p+1)(q+1)} \int \limits_0^\infty \mathrm{L}_p (t) \mathrm{L}_q (t) \mathrm{e}^{-t} \, \mathrm{d} t = \sum \limits_{p=0}^\infty \frac{1}{(p+1)^2} = \zeta(2) = \frac{\pi^2}{6} \, .$$ Del mismo modo, tenemos $$ I = 4\sum \limits_{p=0}^\infty \sum \limits_{q=0}^\infty \sum \limits_{r=0}^\infty \frac{1}{(p+1)(q+1)(r+1)} \int \limits_0^\infty \mathrm{L}_p (t) \mathrm{L}_q (t) \mathrm{L}_r (t) (1- \mathrm{e}^{-t}) \mathrm{e}^{-t} \, \mathrm{d} t \, .$$ Parece que se conocen fórmulas generales para las integrales que implican tres polinomios de Laguerre (véase este papel o este una para una generalización). No sé si son lo suficientemente agradables como para reducir la serie triple a una representación de $\zeta(3)$ sin embargo.
Nota: : Después de hacer algunos cálculos numéricos, ahora sospecho que la serie triple diverge. Esto se debe probablemente al hecho de que la serie original $(\mathrm{B})$ sólo converge en $L^2$ por lo que no se puede utilizar aquí. Sin embargo, para la integral más simple todo funciona.
Por supuesto, es posible combinar los dos métodos aplicando la serie de Laguerre $(\mathrm{B})$ en la ecuación $(\mathrm{A})$ . No sé si estas ideas se pueden utilizar para obtener el resultado final, pero tal vez pueden ayudar a alguien a encontrar una manera.
