Closed-form of $\int_0^1 \operatorname{Li}_3^3(x)\,dx$ and $\int_0^1 \operatorname{Li}_3^4(x)\,dx$

Numerically, $$ I_3 = 27 \zeta (5) \zeta (2)-3 \zeta (3)^2 \zeta (2)+156 \zeta (3) \zeta (2)-\tfrac{153}{8} \zeta (7)-90 \zeta (5)+\zeta (3)^3+21 \zeta (3)^2-660 \zeta (3)-420 \zeta (2)-315 \zeta (4)-\tfrac{477}{4} \zeta (6)+1680 $$ and $I_4$ doesn't seem to have a similar closed form.

The same approach that I used in this answer to decompose the integral as a sum of multiple zeta values works here as well because of the integral form $$ \mathrm{Li}_3(x) = \int_0^x \frac{dt}{t} \int_0^t \frac{du}{u} \int_0^u \frac{du}{1-u}, $$ so it's possible to write $\int_0^1 \mathrm{Li}_3(x)^4\,dx$ as an iterated 13-fold integral, but that looks a little tedious to apply it by hand.

Here's my approach. To reduce clutter let

$$f_m(x) = Li_m(x)$$

and let

$$f_M(x) = f_{m_1}(x)\dots f_{m_n}(x)$$

where $M$ has $n$ entries. Then let

$$I_k(M) = \int f_M(x) dx$$

where we can interpret this as an indefinite or definite integral. Consider the following derivative

$$(x^{k+1}f_M)' = (k+1)x^k f_M + x^{k+1} \sum_{k=1}^n f_{m_k}'\prod_{i\neq k}f_{m_i} = (k+1)x^k f_M + x^{k} \sum_{k=1}^n f_{m_k-1}\prod_{i\neq k}f_{m_i}$$

Integrating with respect to $x$ (and ignoring the extra constant) yields

$$x^{k+1}f_M = (k+1)I_k(M) + \sum_{i=1}^n I_k(M)_i$$

where the notation $I_k(M)_i$ means we subtract $1$ from the $i$-th entry in $M$. This works for every $k+1 \neq 0$. For the case where $k+1 = 0$ we have

$$f_M' = \sum_{k=1}^n f_{m_k}'\prod_{i\neq k}f_{m_i} = \frac{1}{x}\sum_{k=1}^n f_{m_k-1}\prod_{i\neq k}f_{m_i}$$

Integrating yields the formula

$$f_M = \sum_{i=1}^n I_{-1}(M)_i$$

The integrals you seek are $I_0(3,3,3)$ and $I_0(3,3,3,3)$. The difficulty is in the number of terms and not necessarily the size of each term. We can write the recurrence relation

$$I_0(M) = x f_M - \sum_{i=1}^n I_0(M)_i$$

These identities should also hold for the definite version of the integrals.

Anyway, we can start with $I_0(3,3,3)$ and reduce it to $I_0(2,3,3)$. Then we reduce it further to $I_0(1,3,3)$ and $I_0(2,2,3)$ and so on. The question is where do we stop. One should ideally stop when one can compute all the quantities. If not, then one should go further. In the end I think stopping when the smallest entry is 1 might be a good strategy. We reduce

$$(3,3,3) \to (2,3,3) \to (1,3,3),(2,2,3)$$

and then

$$(2,2,3) \to (1,2,3),(2,2,2)$$

The definite version of $I_0(2,2,2)$ corresponds to $\mathcal{I}_3$ from Kirill's linked answer so we have

$$I_0(3,3,3) = 6xf_2^2 f_3- 3xf_2 f_3^2 + x f_3^3 - 6 I_0(2,2,2) - 12 I_0(1,2,3) + 3 I_0(1,3,3)$$

Letting $x \to 1$ and $x \to 0$ we find

$$I_0(3,3,3) = 6 \zeta(2)^2 \zeta(3) - 3 \zeta(2) \zeta(3)^2 + \zeta(3)^3 - 6 \mathcal{I_3} - 12 I_0(1,2,3) + 3 I_0(1,3,3)\\ = 540 - 108 \zeta(2) - 216 \zeta(3) + 72 \zeta(2) \zeta(3) + 6 \zeta(2)^2 \zeta(3) - 3 \zeta(3) \zeta(3)^2 + \zeta(3)^3 - 171 \zeta(4) - 36 \zeta(5) - \tfrac{105}{4} \zeta(6) - 12 I_0(1, 2, 3) + 3 I_0(1, 3, 3)$$

Therefore the problem reduces to calculating The quantities

$$I_0(1, 3, 3) = \int_0^1 Li_1(x)Li_3(x)^2 dx = - \int_0^1 \log(1-x)Li_3(x)^2 dx$$

$$I_0(1, 2, 3) = \int_0^1 Li_1(x)Li_2(x)Li_3(x) dx = -\int_0^1 \log(1-x)Li_2(x)Li_3(x) dx$$

I'm not sure how much of a reduction this is because the number of terms has not been decreased. However, from Kirill's answer, it seems that it'll make his approach easier since one would need fewer iterated integrals.

EDIT: As Kirill mentioned, the procedure reduces the overall weight $|M| = m_1+\dots+m_n$ which could be useful on its own. If we allowed the appearance of $0$'s in $M$ we get divergent integrals, so one would have to handle them in some more special way. If we just went ahead with it we'd get

$$I_0(3,3,3) = 90 f_1^3 x+18 \left(f_3-5 f_2\right) f_1^2 x+3 \left(12 f_2^2-6 f_3f_2+f_3^2\right) f_1 x-6f_2^3 x+f_3^3 x-3 f_2 f_3^2 x+6 f_2^2 f_3 x + R$$

where $R$ consists of

$$R = -270 I_0(0, 1, 1) + 180 I_0(0, 1, 2)- 36 I_0(0, 1, 3) - 36 I_0(0, 2, 2) + 18 I_0(0, 2,3) - 3 I_0(0, 3, 3)$$

We've reduced the weight from 9 to 6, 5, 4, 3 and 2. One could even go a little further and think about making the weight 0, although I'm not sure if that would be a good way to go about this. In any case, we can try this and we obtain a larger sum with $M$'s $(-6,3,3)$, $(-5,2,3)$, $(-4,1,3)$, $(-4,2,2)$, $ (-3,0,3)$, $(-3,1,2)$, $(-2,-1,3)$, $(-2,0,2)$, $(-2,1,1)$, $(-1,-1,2)$, $(-1,0,1)$ and $(0,0,0)$.