Integrals of integer powers of dilogarithm function

I'm interested in evaluating integrals of positive integer powers of the dilogarithm function. I'd like to see the general case tackled if possible, or barring that then as many particular cases as possible.

Problem. For each $n\in\mathbb{N}$, evaluate the definite integral $\mathcal{I}_n$: $$\mathcal{I}_n:=\int_{0}^{1}\left[\operatorname{Li}_2{(x)}\right]^n\,\mathrm{d}x=\,???$$

The first two cases aren't that difficult to evaluate, but for reference I'll just state their values without proof:

$$\begin{align} \int_{0}^{1}\operatorname{Li}_2{(x)}\,\mathrm{d}x &=\zeta{(2)}-1\\ &=0.6449340668482264\dots; \end{align}$$

and

$$\begin{align} \int_{0}^{1}\operatorname{Li}_2{(x)}^2\,\mathrm{d}x &=\frac52\zeta{(4)}-4\zeta{(3)}-2\zeta{(2)}+6\\ &=0.6077123379430154\dots. \end{align}$$

The next case, $n=3$, is much harder than the previous two. However, by combining my answer to this question with Omran Kouba's answer to the same question, and assuming both answers are correct, I was able to infer that

$$\zeta{(5)}+6\zeta{(3)}+\frac{\pi^4}{15}-15=\frac{1}{3!}\int_{0}^{1}\left[\operatorname{Li}_2{(x)}-\zeta{(2)}\right]^3\,\mathrm{d}x.$$

Then by using the Binomial Theorem to expand the integrand on the RHS and by using the previous values of $\mathcal{I}_n$ for $n<3$, I was able to solve for the integral $\mathcal{I}_3$ and find the following value:

$$\begin{align} \mathcal{I}_3 &=\int_{0}^{1}\operatorname{Li}_2{(x)}^3\,\mathrm{d}x\\ &=6\zeta{(5)}+36\zeta{(3)}+\frac{\pi^6}{216}+\frac{19\pi^4}{60}+3\pi^2-2\pi^2\zeta{(3)}-90\\ &=0.6738641012555397\dots. \end{align}$$

But surely there is a less circuitous way to prove this. Can anybody offer a more direct proof that the integral $\mathcal{I}_3$ has the conjectured value indicated above?

Prove: $$\int_{0}^{1}\operatorname{Li}_2{(x)}^3\,\mathrm{d}x\stackrel{?}{=}6\zeta{(5)}+36\zeta{(3)}+\frac{\pi^6}{216}+\frac{19\pi^4}{60}+3\pi^2-2\pi^2\zeta{(3)}-90.$$

And finally, what of the $n>3$ cases? What is the largest $n$ for which $\mathcal{I}_n$ can be evaluated? Can $\mathcal{I}_4$ be evaluated? At the bottom of Omran Kouba's response to the question I referred to above, he says that, to his knowledge, the values of the integral of integer powers of the dilogarithm for powers larger than $2$ are not known. Clearly, for that statement to be true then, at the very least, it would have to be amended to reflect the fact that the value of the integral of the third power of the dilogarithm is indeed known (assuming I'm not somehow the first person to know it!). But are the values for $n>3$ really unknown?

Update: As SuperAbound was able to show, finding $\int_{0}^{1}\operatorname{Li}_2{(x)}^3\,\mathrm{d}x$ via repreated integration by parts actually isn't that difficult at all! For comparison, let's see how much harder finding $\int_{0}^{1}\operatorname{Li}_2{(x)}^4\,\mathrm{d}x$ is if we try to recycle the same method.

To evaluate $\int_0^1\operatorname{Li}_2{(x)}^4\mathrm{d}x$, integrate by parts using

$$\frac{d}{dx}\operatorname{Li}_2{(x)}^{4}=-\frac{4\ln{(1-x)}\operatorname{Li}_2{(x)}^3}{x}.$$

Then,

$$\begin{align} \mathcal{I}_4 &=\int_{0}^{1}\operatorname{Li}_2{(x)}^4\,\mathrm{d}x\\ &=\left[x\operatorname{Li}_2{(x)}^4\right]_{0}^{1}-\int_{0}^{1}x\left(-\frac{4\ln{(1-x)}\operatorname{Li}_2{(x)}^3}{x}\right)\,\mathrm{d}x\\ &=\operatorname{Li}_2{(1)}^4+4\int_{0}^{1}\ln{(1-x)}\operatorname{Li}_2{(x)}^3\,\mathrm{d}x\\ &=\zeta{(2)}^4+4\int_{0}^{1}\ln{(1-x)}\operatorname{Li}_2{(x)}^3\,\mathrm{d}x. \end{align}$$

Next, evaluate $\int_{0}^{1}\ln{(1-x)}\operatorname{Li}_2{(x)}^3\,\mathrm{d}x$ by integrating by parts again using

$$\int\ln{(1-x)}\,\mathrm{d}x=(x-1)\ln{(1-x)}-x+constant,$$

$$\frac{d}{dx}\operatorname{Li}_2{(x)}^3=-\frac{3\ln{(1-x)}\operatorname{Li}_2{(x)}^2}{x}.$$

Then,

$$\begin{align} \int_{0}^{1}\ln{(1-x)}\operatorname{Li}_2{(x)}^3\,\mathrm{d}x &=\left[\left((x-1)\ln{(1-x)}-x\right)\operatorname{Li}_2{(x)}^3\right]_{0}^{1}\\ &~~~~~-\int_{0}^{1}\left((x-1)\ln{(1-x)}-x\right)\left(-\frac{3\ln{(1-x)}\operatorname{Li}_2{(x)}^2}{x}\right)\,\mathrm{d}x\\ &=-\operatorname{Li}_2{(1)}^3\\ &~~~~~-3\int_{0}^{1}\left[\frac{\ln^2{(1-x)}}{x}-\ln^2{(1-x)}+\ln{(1-x)}\right]\operatorname{Li}_2{(x)}^2\,\mathrm{d}x\\ &=-\zeta{(2)}^3-3\int_{0}^{1}\frac{\ln^2{(1-x)}\operatorname{Li}_2{(x)}^2}{x}\,\mathrm{d}x\\ &~~~~~+3\int_{0}^{1}\ln^2{(1-x)}\operatorname{Li}_2{(x)}^2\,\mathrm{d}x-3\int_{0}^{1}\ln{(1-x)}\operatorname{Li}_2{(x)}^2\,\mathrm{d}x. \end{align}$$

The third integral in the last line above has the value,

$$\int_{0}^{1}\ln{(1-x)}\operatorname{Li}_2{(x)}^2\,\mathrm{d}x=2\zeta{(5)}+\frac{19}{2}\zeta{(4)}+12\zeta{(3)}+6\zeta{(2)}-4\zeta{(3)}\zeta{(2)}-30.$$

(to be continued...)

There is a reasonably easy way to do these integrals. Above all, it is systematic. Each dilogarithm can be written as an iterated integral $$ L(x) = \int_0^x \frac{dt}{t}\int_0^t\frac{du}{1-u}. $$ Every power of a dilogarithm can then be expanded into a long product of integrals, giving, as an example: $$ \int_0^1 L(x)^2\,dx = \int_0^1dx \int_0^x \frac{dt_1}{t_1} \int_0^{t_1}\frac{du_1}{1-u_1} \int_0^x \frac{dt_2}{t_2} \int_0^{t_2}\frac{du_2}{1-u_2}. $$ The reason this is easy is that one integral representation of multiple zeta values has the form $$ \int_0^1 \frac{dw_1}{w_1-b_1}\int_0^{w_1}\frac{dw_2}{w_2-b_2}\cdots\int_0^{w_{n-1}}\frac{dw_n}{w_n-b_n}, $$ where $b_j\in\{0,1\}$.

In the case of $\int L(x)^n$, the domain of the integral is $$ 0< t_1,\ldots,t_n < x < 1, \qquad 0 < u_j < t_j, $$ which can be rearranged into a union of simplices of the form $$ 0<w_{2n}<\cdots<w_{1}<x<1, $$ (where $w_j$'s are all of $t_j$'s and $u_j$'s in some order), and the integral over each simplex has the form $$ \int_0^1dx\,\Omega, $$ where $\Omega$ is an iterated integral expression, starting with either $\int_0^x \frac{dw_1}{w_1}$ or $\int_0^x \frac{dw_1}{1-w_1}$.

It is easy to interchange the order of integration of $x$ and $w_1$, integrate over $x$, and be left with a very similar integral over another simplex with only the variables $w_2,\ldots,w_{2n}$ present.

Implementing this recursive procedure gives the following expressions for $I$: $$\begin{eqnarray} I_1 &=& \zeta(2) - 1 \\I_2 &=& -2\zeta(2) - 4\zeta(2,1) + 2\zeta(2,2) + 4\zeta(3,1) + 6 \\I_3 &=& 18\zeta(2) + 36\zeta(2,1) + 36\zeta(2,1,1) - 12\zeta(2,1,2) - 6\zeta(2,2) - 24\zeta(2,2,1) + 6\zeta(2,2,2) + 12\zeta(2,3,1) - 12\zeta(3,1) - 36\zeta(3,1,1) + 12\zeta(3,1,2) + 24\zeta(3,2,1) + 36\zeta(4,1,1) - 90 \\I_4 &=& -360\zeta(2) - 720\zeta(2,1) - 864\zeta(2,1,1) - 576\zeta(2,1,1,1) + 144\zeta(2,1,1,2) + 144\zeta(2,1,2) + 288\zeta(2,1,2,1) - 48\zeta(2,1,2,2) - 96\zeta(2,1,3,1) + 72\zeta(2,2) + 288\zeta(2,2,1) + 432\zeta(2,2,1,1) - 96\zeta(2,2,1,2) - 24\zeta(2,2,2) - 192\zeta(2,2,2,1) + 24\zeta(2,2,2,2) + 48\zeta(2,2,3,1) - 48\zeta(2,3,1) - 288\zeta(2,3,1,1) + 48\zeta(2,3,1,2) + 96\zeta(2,3,2,1) + 144\zeta(2,4,1,1) + 144\zeta(3,1) + 432\zeta(3,1,1) + 576\zeta(3,1,1,1) - 144\zeta(3,1,1,2) - 48\zeta(3,1,2) - 288\zeta(3,1,2,1) + 48\zeta(3,1,2,2) + 96\zeta(3,1,3,1) - 96\zeta(3,2,1) - 432\zeta(3,2,1,1) + 96\zeta(3,2,1,2) + 192\zeta(3,2,2,1) + 288\zeta(3,3,1,1) - 144\zeta(4,1,1) - 576\zeta(4,1,1,1) + 144\zeta(4,1,1,2) + 288\zeta(4,1,2,1) + 432\zeta(4,2,1,1) + 576\zeta(5,1,1,1) + 2520 \end{eqnarray}$$

I computed these values up to $n=8$, but the expressions are long ($I_5$ has 428 terms), and the integers involved become prohibitively large, so I didn't include them here.

For small $n$, at least up to $n=4$, I found numerically that the above expressions can be written in terms of ordinary zeta values, as follows: $$\begin{eqnarray}\def\tfrac#1#2{{\textstyle\frac{#1}{#2}}} I_1 &=& \zeta(2) - 1 \\I_2 &=& -4 \zeta (3)-2 \zeta (2)+\tfrac{5}{2} \zeta (4)+6 \\I_3 &=& -12 \zeta (3) \zeta (2)+6 \zeta (5)+36 \zeta (3)+18 \zeta (2)+\tfrac{57}{2} \zeta (4)+\tfrac{35}{8} \zeta (6)-90 \\I_4 &=& 24 \zeta (5) \zeta (2)+144 \zeta (3) \zeta (2)+84 \zeta (3) \zeta (4)-195 \zeta (7)-648 \zeta (5)+144 \zeta (3)^2-720 \zeta (3)-360 \zeta (2)-774 \zeta (4)-\frac{11}{2} \zeta (6)+\frac{175}{24} \zeta (8)+2520 \end{eqnarray}$$ I checked $n\leq3$ using known identities satisfied by MZVs, such as those here. For $I_4$, while this is not a proof that the expression is correct, this entire procedure reduces the problem of evaluating the integrals $I_n$ to the problem of applying known identities to multiple zeta values, which is more straightforward.

A quick integer relation search for a closed form for $I_5$ gave me nothing.

There may or may not be a closed form for $I_5$ (I believe this is not a settled question), but there is no systematic way to reduce multiple zeta values to ordinary zeta values. For example, weight 10, which is what $I_5$ requires, is the first weight at which known relations between depth-two zeta values are insufficient to reduce them to ordinary zeta values ("Euler Sums and Contour Integral Representations" by Flajolet and Salvy). Because here the depth can be up to five, this problem is harder.

Update (by editor after 6 years): Settled now. By further computation one have

$$\small\int_0^1 \text{Li}_2(x){}^5 \, dx=1440 \zeta(6,2)+240 \zeta (3)^3+120 \pi ^2 \zeta (3)^2-4320 \zeta (3)^2-\frac{23 \pi ^6 \zeta (3)}{189}-110 \pi ^4 \zeta (3)-600 \pi ^2 \zeta (3)-9360 \zeta (5) \zeta (3)+25200 \zeta (3)-\frac{127 \pi ^4 \zeta (5)}{3}-540 \pi ^2 \zeta (5)+30600 \zeta (5)-\frac{325 \pi ^2 \zeta (7)}{2}+14445 \zeta (7)+\frac{16495 \zeta (9)}{3}+\frac{\pi ^{10}}{7776}+\frac{51137 \pi ^8}{45360}+\frac{1907 \pi ^6}{126}+335 \pi ^4+2100 \pi ^2-113400$$

Here $ \zeta(6,2)$ is irreducible. In fact all $I_n$ can be reduced to MZVs. The generalized algorithm I used also gives the proof to closed-form of $I_4$ above.

Allow me to present another approach. Though lengthier and more tedious, this method is independent of results derived from previous answers. The main idea is to expand ${\rm Li}_2(x)$ as a series and integrate term by term.

For $n=1$, \begin{align} \mathcal{I}_1 &=\int^1_0{\rm Li}_2(x) \ {\rm d}x\\ &=\int^1_0\sum^\infty_{j=1}\frac{x^j}{j^2}{\rm d}x\\ &=\sum^\infty_{j=1}\frac{1}{j^2(j+1)}\\ &=\sum^\infty_{j=1}\frac{1}{j^2}-\frac{1}{j}+\frac{1}{j+1}\\ \end{align} We see that \begin{align} \sum^\infty_{k=1}\frac{1}{k^2} &=\frac{\pi^2}{6}\\ \sum^\infty_{k=1}\left(\frac{1}{k+1}-\frac{1}{k}\right) &=\sum^\infty_{k=2}\frac{1}{k}-\sum^\infty_{k=1}\frac{1}{k}\\ &=-1 \end{align} Therefore $$\mathcal{I}_1=\frac{\pi^2}{6}-1$$

For $n=2$, \begin{align} \mathcal{I}_2 &=\int^1_0{\rm Li}^2_2(x) \ {\rm d}x\\ &=\int^1_0\sum^\infty_{j=1}\frac{x^j}{j^2}\sum^\infty_{k=1}\frac{x^k}{k^2}{\rm d}x\\ &=\sum^\infty_{j=1}\frac{1}{j^2}\sum^\infty_{k=1}\frac{1}{k^2(k+j+1)}\\ &=\sum^\infty_{j=1}\frac{1}{j^2(j+1)}\sum^\infty_{k=1}\frac{1}{k^2}-\frac{1}{k(j+1)}+\frac{1}{(k+j+1)(j+1)}\\ &=\frac{\pi^2}{6}\left(\frac{\pi^2}{6}-1\right)+\sum^\infty_{j=1}\frac{1}{j^2(j+1)^2}\sum^\infty_{k=1}\frac{1}{k+j+1}-\frac{1}{k}\\ &=\frac{\pi^4}{36}-\frac{\pi^2}{6}+\color{red}{\sum^\infty_{j=1}\frac{2H_{j+1}}{j}}\color{green}{-\sum^\infty_{j=1}\frac{2H_{j+1}}{j+1}}\color{blue}{-\sum^\infty_{j=1}\frac{H_{j+1}}{j^2}}\color{purple}{-\sum^\infty_{j=1}\frac{H_{j+1}}{(j+1)^2}}\\ \end{align} It is trivial to show that \begin{align} \sum^\infty_{j=1}\frac{H_j}{j^2} &=2\zeta(3) \end{align} We see that \begin{align} \color{red}{\sum^\infty_{j=1}\frac{2H_{j+1}}{j}} &=\sum^\infty_{j=1}\frac{2H_j}{j}+\frac{2}{j(j+1)}\\ &=\color{red}{\sum^\infty_{j=1}\frac{2H_j}{j}+2}\\ \color{green}{-\sum^\infty_{j=1}\frac{2H_{j+1}}{j+1}} &=\color{green}{-\sum^\infty_{j=1}\frac{2H_j}{j}+2}\\ \color{blue}{-\sum^\infty_{j=1}\frac{H_{j+1}}{j^2}} &=-\sum^\infty_{j=1}\frac{H_{j}}{j^2}-\sum^\infty_{j=1}\frac{1}{j^2(j+1)}\\ &=\color{blue}{-2\zeta(3)+1-\frac{\pi^2}{6}}\\ -\color{purple}{\sum^\infty_{j=1}\frac{H_{j+1}}{(j+1)^2}} &=-\sum^\infty_{j=1}\frac{H_{j}}{j^2}+1\\ &=\color{\purple}{-2\zeta(3)+1} \end{align} Therefore $$\mathcal{I}_2=\frac{\pi^4}{36}-4\zeta(3)-\frac{\pi^2}{3}+6$$

For $n=3$, \begin{align} \mathcal{I}_3 &=\int^1_0{\rm Li}^3_2(x) \ {\rm d}x\\ &=\int^1_0\sum^\infty_{j=1}\frac{x^j}{j^2}\sum^\infty_{k=1}\frac{x^k}{k^2}\sum^\infty_{m=1}\frac{x^m}{m^2}{\rm d}x\\ &=\sum^\infty_{j=1}\frac{1}{j^2}\sum^\infty_{k=1}\frac{1}{k^2}\sum^\infty_{m=1}\frac{1}{m^2(m+j+k+1)}\\ &=\sum^\infty_{j=1}\frac{1}{j^2}\sum^\infty_{k=1}\frac{1}{k^2(j+k+1)}\sum^\infty_{m=1}\frac{1}{m^2}-\frac{1}{m(j+k+1)}+\frac{1}{(m+j+k+1)(j+k+1)}\\ &=\frac{\pi^2}{6}\mathcal{I}_2-\sum^\infty_{j=1}\frac{1}{j^2}\sum^\infty_{k=1}\frac{H_{k+j+1}}{k^2(k+j+1)^2}\\ &=\frac{\pi^2}{6}\mathcal{I}_2-\sum^\infty_{j=1}\frac{1}{j^2(j+1)^2}\sum^\infty_{k=1}\left(\frac{H_{k+j+1}}{k^2}+\frac{H_{k+j+1}}{(k+j+1)^2}\right)+2\sum^\infty_{j=1}\frac{1}{j^2(j+1)^3}\sum^\infty_{k=1}\left(\frac{H_{k+j+1}}{k}-\frac{H_{k+j+1}}{k+j+1}\right) \end{align} The first inner sum is \begin{align} \sum^\infty_{k=1}\frac{H_{k+j+1}}{k^2} &=2\zeta(3)+\sum^\infty_{k=1}\sum^{j+1}_{m=1}\frac{1}{k^2(k+m)}\\ &=2\zeta(3)+\sum^\infty_{k=1}\sum^{j+1}_{m=1}\frac{1}{mk^2}-\frac{1}{m^2k}+\frac{1}{m^2(k+m)}\\ &=2\zeta(3)+\frac{\pi^2}{6}H_{j+1}-\sum^{j+1}_{m=1}\frac{1}{m^2}\sum^\infty_{k=1}\left(\frac{1}{k}-\frac{1}{k+m}\right)\\ &=2\zeta(3)+\frac{\pi^2}{6}H_{j+1}-\sum^{j+1}_{k=1}\frac{H_k}{k^2}\\ \end{align} The second sum is \begin{align} \sum^\infty_{k=1}\frac{H_{k+j+1}}{(k+j+1)^2} &=2\zeta(3)-\sum^{j+1}_{k=1}\frac{H_k}{k^2} \end{align} Adding them up gives $$4\zeta(3)+\frac{\pi^2}{6}H_{j+1}-2\sum^{j+1}_{k=1}\frac{H_k}{k^2}$$ The third sum is \begin{align} \sum^\infty_{k=1}\frac{H_{k+j+1}}{k} &=\sum^\infty_{k=1}\frac{H_k}{k}+\sum^\infty_{k=1}\sum^{j+1}_{m=1}\frac{1}{k(k+m)}\\ &=\sum^\infty_{k=1}\frac{H_{k}}{k}+\sum^\infty_{k=1}\sum^{j+1}_{m=1}\frac{1}{mk}-\frac{1}{m(k+m)}\\ &=\sum^\infty_{k=1}\frac{H_{k}}{k}+\sum^{j+1}_{k=1}\frac{H_k}{k}\\ \end{align} The fourth sum is \begin{align} \sum^\infty_{k=1}\frac{H_{k+j+1}}{k+j+1} &=\sum^\infty_{k=1}\frac{H_{k}}{k}-\sum^{j+1}_{k=1}\frac{H_k}{k} \end{align} Subtracting them gives $$H_{j+1}^2+H_{j+1}^{(2)}$$ So $$\mathcal{I}_3=\frac{\pi^2}{6}\mathcal{I}_2-\sum^\infty_{j=1}\frac{1}{j^2(j+1)^2}\left(4\zeta(3)+\frac{\pi^2}{6}H_{j+1}-2\sum^{j+1}_{k=1}\frac{H_k}{k^2}\right)+2\sum^\infty_{k=1}\frac{H_{j+1}^2+H_{j+1}^{(2)}}{j^2(j+1)^3}$$ I will try to complete the evaluation of $\mathcal{I}_3$ over the next few days, and add in the computation of $\mathcal{I}_4$ using this method if I have the time. It is clear that $\mathcal{I}_3$ can be reduced to a bunch of Euler sums, which are all rather manageable to compute.

Integrate by parts twice. \begin{align} \int^1_0{\rm Li}_2^3(x){\rm d}x &=\left[x{\rm Li}^3_2(x)\right]^1_0+\int^1_03{\rm Li}^2_2(x)\ln(1-x){\rm d}x\\ &=\frac{\pi^6}{216}-\frac{\pi^4}{12}+6\color\red{\int^1_0\frac{\left[(x-1)\ln(1-x)-x\right]{\rm Li}_2(x)\ln(1-x)}{x}{\rm d}x}\\ &=\frac{\pi^6}{216}-\frac{\pi^4}{12}+6\zeta(5)+\frac{2\pi^4}{5}+36\zeta(3)+3\pi^2-2\pi^2\zeta(3)-90\\ &=\frac{\pi^6}{216}+\frac{19\pi^4}{60}+6\zeta(5)+36\zeta(3)+3\pi^2-2\pi^2\zeta(3)-90\\ \end{align} The red integral was evaluated in your previous answer.