Conjectured closed form for $\int_0^1\frac{\operatorname{li}^4(x)}{x^4}\,dx$
Solution 1:
This is not a complete answer, but it gives an alternative representation of $I(4)$ which may be better to work with. Define $F(s)$ by $$ F(s) = \int_{0}^{1} \frac{s}{1-sx} \log(1-x)\operatorname{Li}_{2}(x) \, dx. $$
Then we have the following representation:
Claim. We have $\displaystyle I(4) = \frac{71\pi^4}{540} - \frac{8}{3}F(1/2)$.
Consequently, the conjecture is equivalent to proving
$$ \color{red}{ F(1/2) \stackrel{?}{=} 4\operatorname{Li}_4 (1/2) + \frac{1}{6}\log^4 2 - \zeta(2)\log^2 2 - \frac{11}{480} \pi^4. } $$
These coefficients are sort of familiar, so I believe that we are on a right track. On the other hand, this still seems daunting because
$$ F(1/2) = \int_{0}^{1} \frac{\log x}{1+x} \left( \zeta(2) - \log x \log(1-x) - \operatorname{Li}_2(x) \right) \, dx $$
is obviously related to the class of alternating Euler sums.
Step 1. Using @Kirill's calculation, we obtain
$$ I(n) = \int_{[0,\infty)^n} \frac{dt_1 \cdots dt_n}{(1+t_1)\cdots(1+t_n)(1+t_1+\cdots+t_n)}. $$
Step 2. We consider the following change of variables:
$$ (x, y) = \Phi(u, v) = (\tfrac{u+v}{2}, \sqrt{uv}). $$
Then $\Phi$ maps the region $\Delta = \{ (u, v) : 0 \leq v \leq u \}$ onto $\Delta$ itself, and its Jacobian is
$$ \frac{\partial(x,y)}{\partial(u,v)} = \frac{2y}{\sqrt{x^2 - y^2}}. $$
From the symmetry and this change of variables, we have
\begin{align*} I(2n) &= \int\limits_{(u_1, v_1), \cdots, (u_n, v_n) \in \Delta} \Bigg[ \prod_{i=1}^{n} \frac{2 du_i dv_i}{(1+u_i)(1+v_i)} \Bigg] \frac{1}{1+u_1+v_1+\cdots+u_n+v_n} \\ &= \int\limits_{(x_1, y_1), \cdots, (x_n, y_n) \in \Delta} \Bigg[ \prod_{i=1}^{n} \frac{4y_i \, dx_i dy_i}{(1+2x_i+y_i^2)\sqrt{x_i^2 - y_i^2}} \Bigg] \frac{1}{1+2x_1 + \cdots + 2x_n}. \end{align*}
We remark that this trick is essentially the same as the method used by @Random Variable. Now using the following formula
$$ \int_{0}^{x} \frac{2y \, dy}{(1+2x+y^2)\sqrt{x^2-y^2}} = \frac{\log(1+2x)}{1+x}, $$
we can make a further simplification and we get
$$ I(2n) = 2^n \int_{[0,\infty)^n} \Bigg[ \frac{\log(1+2x_i) \, dx_i}{1+x_i} \Bigg] \frac{1}{1+2x_1 + \cdots + 2x_n}. $$
In particular, when $n = 2$, then the substitutions $1+2x_1 = \frac{1}{1-x}$ and $1+2x_2 = \frac{1}{1-y}$ yield
$$ \color{blue}{ I(4) = \int_{0}^{1}\int_{0}^{1} \frac{\log(1-x)\log(1-y)}{(1-xy)(1-\frac{1}{2}x)(1-\frac{1}{2}y)} \, dxdy. } \tag{1} $$
This is a reminiscence of the integral
$$ \int_{0}^{1} \int_{0}^{1} \frac{\log(1-x)\log(1-y)}{1-xy} \, dxdy = \frac{17\pi^4}{360}, $$
which already required some hard works as in here. Indeed, this result will be used in our computation.
Step 3. We first remark the following identity:
Lemma 1. For any $s \in \Bbb{C}\setminus[1,\infty)$ and $m = 1, 2, \cdots$, we have $$ \int_{0}^{1} \frac{s}{1-sx} \log^{m-1}(1-x) \, dx = (-1)^{m}\Gamma(m) \operatorname{Li}_{m}\left( \frac{s}{s-1} \right). $$
(I acknowledge that this lemma is motivated by @Felix Marin's calculation.) This is easily proved by letting $w = s/(s-1)$ and noticing that
$$ \frac{s}{1-s x} = - \frac{w}{1-w(1-x)}. $$
To accompany this result, we remind some identities involving dilogarithm:
Lemma 2. For $s \in \Bbb{C} \setminus [1, \infty)$ we have $\displaystyle \operatorname{Li}_{2}\left( \frac{s}{s-1} \right) = -\operatorname{Li}_2(s) - \frac{1}{2}\log^2(1-s)$.
Then from Lemma 1 and Lemma 2, it is straightforward to check that \begin{align*} &\int_{0}^{1} \int_{0}^{1} \frac{\log(1-x)\log(1-y)}{(1-xy)(1-sx)} \, dxdy \\ &= \int_{0}^{1} \frac{\log(1-x)}{x(1-sx)}\operatorname{Li}_{2}\left( \frac{x}{x-1} \right) \, dx \\ &= - \int_{0}^{1} \left( \frac{1}{x} + \frac{s}{1-sx} \right)\left( \operatorname{Li}_{2}(x) + \frac{1}{2}\log^{2}(1-x) \right) \log(1-x) \, dx \\ &= \frac{17\pi^{4}}{360} - F(s) - 3 \operatorname{Li}_{4}\left( \frac{s}{s-1} \right). \end{align*}
Finally, in order to facilitate this calculation, we notice that
\begin{align*} &\frac{1}{(1-xy)(1-sx)(1-sy)} \\ &\qquad = \frac{1}{1-s^{2}}\left( \frac{1}{1-xy}\left( \frac{1}{1-sx} + \frac{1}{1-sy} - 1 \right) - \frac{s}{1-sx} \frac{s}{1-sy} \right). \end{align*}
Consequently it gives
\begin{align*} &\int_{0}^{1} \int_{0}^{1} \frac{\log(1-x)\log(1-y)}{(1-xy)(1-sx)(1-sy)} \, dxdy \\ &= \color{blue}{ \frac{1}{1-s^{2}}\left( \frac{17\pi^{4}}{360} - 2F(s) - 6 \operatorname{Li}_{4}\left( \frac{s}{s-1} \right) - \operatorname{Li}_{2}\left( \frac{s}{s-1} \right)^{2} \right) }. \tag{2} \end{align*}
Then plugging $s = 1/2$ proves the claim. ////
Solution 2:
Although I'm 7 years late, here is an alternative direct method:
First, we shall introduce a Laplace transform that will assist us:
$$e^{2x} \operatorname{Ei}^2 (-x) = \int_{0}^{\infty} \frac{2e^{-x t}\ln (1+t)}{2+t} \, dt$$
This can be proven by using multiplying the following representation with itself: $$-e^x \operatorname{Ei} (-x) = \int_{0}^{\infty} \frac{e^{-x t}}{t+1} \, dt$$
In a similar fashion
$$e^{4 x} \operatorname{Ei}^4 (-x) = \int_{0}^{\infty} \int_{0}^{\infty} \frac{4 e^{-x(t+u)} \ln (1+t) \ln (1+u)}{(2+t)(2+u)}\, du \, dt$$
$$\begin{align} \implies I = \int_{0}^{\infty} e^{3x} \operatorname{Ei}^4 (-x) \, dx &= \int_{0}^{\infty} \int_{0}^{\infty} \int_{0}^{\infty} \frac{4 e^{-x(t+u+1)} \ln (1+t) \ln (1+u)}{(2+t)(2+u)} \, dx \, du \, dt \\ &= \int_{0}^{\infty} \frac{4 \ln (1+t)}{2+t} \int_{0}^{\infty} \frac{\ln (1+u)}{(2+u)(1+t+u)} \, du \, dt \\ &= \int_{0}^{\infty} \frac{-4\ln (1+t)}{(2+t)(t-1)} \left(\frac{\pi^2}{12} + \operatorname{Li}_2 (-t) \right) \, dt\end{align}$$
Enforce the substitution $t = e^x - 1 \> \implies dt = e^x \, dx$: $$I = -4 \int_{0}^{\infty} \frac{x e^x}{(e^x+1)(e^x-2)}\left(\frac{\pi^2}{12}+\operatorname{Li}_2 (1-e^x)\right) \, dx$$
Define $$A=-\frac{4}{3} \int_{0}^{\infty} \frac{x}{e^x+1} \left(\frac{\pi^2}{12}+\operatorname{Li}_2 (1-e^x)\right) \, dx$$ $$B=-\frac{8}{3} \int_{0}^{\infty} \frac{x}{e^x-2} \left(\frac{\pi^2}{12}+\operatorname{Li}_2 (1-e^x)\right) \, dx$$ such that $I=A+B$
Using $$(1-2^{1-s}) \zeta (s) \Gamma (s)=\int_{0}^{\infty} \frac{x^{s-1}}{e^x+1} \, dx$$ $$\operatorname{Li}_2 (x) = \int_{0}^{\infty} \frac{x t}{e^t-x} \, dt$$
We have $$\implies A=-\frac{\pi^4}{108} - \frac{4}{3} \int_{0}^{\infty} \frac{x \operatorname{Li}_2 (1-e^x)}{e^x+1} \, dx$$
$$\implies A = -\frac{\pi^4}{108} - \frac{4}{3} \int_{0}^{\infty} \frac{x \operatorname{Li}_2 (1-e^x)}{1-e^x} \, dx - \frac{8}{3} \int_{0}^{\infty} \frac{x}{e^x-2} \left(\frac{\pi^2}{12}+\operatorname{Li}_2 (1-e^x)\right) \, dx \\= -\frac{\pi^4}{108} - \frac{4}{3} \underbrace{\int_{0}^{\infty} \frac{x \operatorname{Li}_2 (1-e^x)}{1-e^x} \, dx}_J + B$$
Using dilogarithm identities $$J = \int_{0}^{\infty} \frac{x \operatorname{Li}_2 (1-e^x)}{1-e^x} \, dx = \int_{0}^{\infty} \frac{x (\operatorname{Li}_2 (1-e^{-x})+ x^2/2)}{e^x-1} \, dx$$ Using $$\zeta(s) \Gamma (s) = \int_{0}^{\infty} \frac{x^{s-1}}{e^x - 1} \, dx$$
$$\implies J = \frac{\pi^4}{30} + \int_{0}^{\infty} x \left(\frac{\pi^2}{6} - \operatorname{Li}_2 (1-e^{-x})\right) \, dx + \int_{0}^{\infty} x^2 \ln (1-e^{-x}) + \frac{x}{2} \ln^2 (1-e^{-x})\, dx$$ With a few applications of integration by parts we conclude $$J = \frac{17 \pi^4}{360} \implies A = -\frac{13 \pi^4}{180}+B$$
$$I = \frac{29\pi^4}{540} - \frac{8}{3} \int_{0}^{\infty} \frac{x\operatorname{Li}_2 (1-e^x)}{e^x+1} \, dx$$ Applying integration by parts we have $$\begin{align} I &= \frac{29\pi^4}{540} - \frac{8}{3} \int_{0}^{\infty}\frac{x}{e^{-x}-1}\left(\operatorname{Li}_2(-e^{-x})-x\ln(1+e^{-x})\right) \, dx \\ &=\frac{29\pi^4}{540}-\frac{8}{3}\left(\frac{\pi^4}{144}-\frac{7}{2}\ln(2)\zeta(3)\right)-\frac{8}{3}\int_{0}^{\infty}\frac{x \operatorname{Li}_2(-e^{-x})}{1-e^{-x}} \, dx\\ &= \frac{19 \pi^4}{540} + \frac{28}{3} \ln (2) \zeta (3)-\frac{8}{3}\underbrace{\int_{0}^{\infty}\frac{x \operatorname{Li}_2(-e^{-x})}{1-e^{-x}} \, dx}_K \end{align}$$
Once again using $$\operatorname{Li}_2 (x) = \int_{0}^{\infty} \frac{x t}{e^t-x} \, dt$$
We have $$K = \int_{0}^{\infty}\frac{x \operatorname{Li}_2(-e^{-x})}{1-e^{-x}} \, dx=-\frac{\pi^2}{3}\int_{0}^{\infty}\frac{t^2}{1+e^t} \, dt+\int_{0}^{\infty}\frac{t^2\ln(1+e^{-t})}{1+e^t} \, dt+\frac{1}{2}\int_{0}^{\infty}\frac{t\ln^2(1+e^{-t})}{1+e^t}\, dt+\int_{0}^{\infty}\frac{t\operatorname{Li}_2\left(\frac{1}{1+e^{-t}}\right)}{1+e^t} \, dt$$
It is not too much work to determine
$$\implies K =-\frac{7\pi^4}{120}-\frac{\pi^2}{8}\ln^2(2)+\frac{\ln^4(2)}{8}+3\operatorname{Li}_4 \left(\frac{1}{2}\right)+\frac{21}{8}\ln(2)\zeta(3)-\int_{0}^{\infty}x\left(\frac{\pi^2}{12}+\operatorname{Li}_2\left(\frac{1}{e^{-x}-1}\right)\right)\, dx$$
$$-\int_{0}^{\infty}x\left(\frac{\pi^2}{12}+\operatorname{Li}_2\left(\frac{1}{e^{-x}-1}\right)\right)\, dx = \int_{0}^{\infty}\frac{x^2}{2(e^x-1)}\left(\ln\left(1-\frac{e^{-x}}{2}\right)-\ln(1-e^{-x})+\ln(2)\right)\, dx$$
I will state a proposition that we shall now use that I will omit the proof of as it is somewhat tedious:
Proposition 1 $$f(b) = \int_{0}^{\infty}\frac{x^2\ln(1+be^{-x})}{2(e^x-1)}dx=-\frac{1}{2}\text{Li}_2(-b)^2-\ln(1+b)\text{Li}_3(-b)+\text{Li}_4(-b)+\ln(1+b)\zeta(3)$$
Putting everything together we have then
$$\implies I = \frac{19 \pi^4}{540} + \frac{28}{3} \ln (2) \zeta (3) -\frac{8}{3}\left(-\frac{7\pi^4}{120}-\frac{\pi^2}{8}\ln^2(2)+\frac{\ln^4(2)}{8}+3\operatorname{Li}_4\left(\frac{1}{2}\right)+\frac{21}{8}\ln(2)\zeta(3)+\ln(2)\zeta(3)+f\left(-\frac{1}{2}\right)-f(-1)\right)$$
Note that $f(-1)$ should be realised as $\lim_{b \, \to \, -1} f(b)$.
Upon simplification, this gives the desired result $$I=\frac{26\pi^4}{135}+\frac49\left(\pi^2\ln^22-\ln^42\right)-\frac{32}{3}\operatorname{Li}_4\!\left(\tfrac12\right)$$ $\square$