Major Edit: This is almost a completely new response rather than an edit. The previous version of this response was extremely long and clumsy, and ultimately failed to even yield a definite final result. This new and improved response is much more streamlined and does include a definite final value.


Evaluation of integral $I_3$:

The triple integral defining $I_3$ can in principle be integrated in any order, but integrating in "alphabetical order" (i.e., with the integral over $x$ as the outermost one, and the integral over $z$ as the innermost) is probably the best way to go and is the order used in the first step below. Next, we rescale the integral over $z$ via the substitution $t=(xy)\,z$; after that, we also rescale the integral over $y$ via the substitution $u=(x)\,y$. Now, instead of evaluating the integrals from innermost-to-outermost, note that our integral is now in a form that lends itself very well to integration by parts with respect to $x$. The result is a sum of two double integrals:

$$\begin{align} I_3 &=\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}\ln{(1-x)}\ln{(1-xy)}\ln{(1-xyz)}\,\mathrm{d}x\mathrm{d}y\mathrm{d}z\\ &=\int_{0}^{1}\mathrm{d}x\,\ln{(1-x)}\int_{0}^{1}\mathrm{d}y\,\ln{(1-xy)} \int_{0}^{1}\mathrm{d}z\,\ln{(1-xyz)}\\ &=\int_{0}^{1}\mathrm{d}x \frac{\ln{(1-x)}}{x} \int_{0}^{x}\mathrm{d}u \frac{\ln{(1-u)}}{u}\int_{0}^{u}\mathrm{d}t\,\ln{(1-t)}\\ &=-\operatorname{Li}_2{(1)}\int_{0}^{1}\mathrm{d}u \frac{\ln{(1-u)}}{u}\int_{0}^{u}\mathrm{d}t\,\ln{(1-t)}\\ &~~~~~+\int_{0}^{1}\mathrm{d}x\frac{\operatorname{Li}_2{(x)}\ln{(1-x)}}{x}\int_{0}^{x}\mathrm{d}t\,\ln{(1-t)}\\ &=:J+K. \end{align}$$

Using the evaluations of $J$ and $K$ below, we arrive at a final value for $I_3$:

$$\begin{align} I_3 &=J+K\\ &=\left[-3\zeta{(2)}+2\zeta{(3)}\zeta{(2)}\right]+\left[\zeta{(5)}+6\zeta{(4)}+6\zeta{(3)}+3\zeta{(2)}-2\zeta{(3)}\zeta{(2)}-15\right]\\ &=\zeta{(5)}+6\zeta{(4)}+6\zeta{(3)}-15\\ &=\zeta{(5)}+6\zeta{(3)}+\frac{\pi^4}{15}-15\\ &=-0.2567 9142 3632 2352\dots . \end{align}$$

$$I_3=\color{blue}{\zeta{(5)}+6\zeta{(3)}+\frac{\pi^4}{15}-15}.$$


Evaluation of integral $J$:

First we state without proof the following three anti-derivatives:

$$\int\mathrm{d}u\ln{(1-u)}=(u-1)\ln{(1-u)}-u+constant;$$

$$\int\mathrm{d}u\ln^2{(1-u)}=(u-1)\left(\ln^2{(1-u)}-2\ln{(1-u)}+2\right)+constant;$$

$$\int\mathrm{d}u\frac{\ln^2{(1-u)}}{u}=-2\operatorname{Li}_3{(1-u)}+2\operatorname{Li}_2{(1-u)}\ln{(1-u)}+\ln{(u)}\ln^2{(1-u)}+constant.$$

They may each be easily verified by differentiating the right-hand-sides, or checked using WolframAlpha. Once obtained, the integral $J$ may be calculated directly from these three anti-derivatives:

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


Evaluation of integral $K$:

$$\begin{align} K &=\int_{0}^{1}\mathrm{d}x\frac{\operatorname{Li}_2{(x)}\ln{(1-x)}}{x}\int_{0}^{x}\mathrm{d}t\,\ln{(1-t)}\\ &=\int_{0}^{1}\mathrm{d}x\frac{\operatorname{Li}_2{(x)}\ln{(1-x)}}{x} \left[(x-1)\ln{(1-x)}-x\right]\\ &=\int_{0}^{1}\mathrm{d}x \operatorname{Li}_2{(x)}\ln{(1-x)} \left[-1+\ln{(1-x)}-\frac{\ln{(1-x)}}{x}\right]\\ &=\int_{0}^{1}\mathrm{d}x \operatorname{Li}_2{(x)} \left[-\ln{(1-x)}+\ln^2{(1-x)}-\frac{\ln^2{(1-x)}}{x}\right]\\ &=-\int_{0}^{1}\mathrm{d}x \ln{(1-x)}\operatorname{Li}_2{(x)} + \int_{0}^{1}\mathrm{d}x \ln^2{(1-x)}\operatorname{Li}_2{(x)} - \int_{0}^{1}\mathrm{d}x \frac{\ln^2{(1-x)} \operatorname{Li}_2{(x)}}{x}\\ &=K_1+K_2+K_3\\ &=\left[2\zeta{(3)}+\zeta{(2)}-3\right]+\left[6\zeta{(4)}+4\zeta{(3)}+2\zeta{(2)}-12\right]+\left[\zeta{(5)}-2\zeta{(3)}\zeta{(2)}\right]\\ &=\zeta{(5)}+6\zeta{(4)}+6\zeta{(3)}+3\zeta{(2)}-2\zeta{(3)}\zeta{(2)}-15. \end{align}$$

The evaluations of $K_1$ and $K_2$ can easily be found by first using a CAS to find the anti-derivatives. Finally, $K_3$ is calculated below.


Evaluation of integral $K_3$:

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

The first integral in the last line above has already been calculated as part of the evaluation of integral $J$. The evaluations of the second and third integrals can be found in the responses to this question and this question, respectively.

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


This is not a solution, but it explains why $I_n$ for $n\geq 3$ is difficult. Indeed, the change of variables $(y_1,\ldots,y_n)=(x_1,x_1x_2,\ldots,x_1x_2\ldots x_n)$, (i.e. $y_k=x_1x_2\ldots x_k$) shows that $$\eqalign{ I_n&=\int_{1\geq y_1\geq y_2\geq \cdots\geq y_n\geq 0}\ln(1-y_1) \ln(1-y_2)\cdots\ln(1-y_n)\frac{dy_1\cdots d y_n}{y_1\cdots y_{n-1}}\cr &=\int_{\color{red}{y_n}=0}^1\left(\int_{1\geq y_1\geq y_2\geq \cdots\geq y_{n-1}\geq \color{red}{y_n}} \prod_{k=1}^{n-1}\frac{\ln(1-y_k)}{y_k}dy_1\ldots dy_{n-1}\right)\ln(1-y_n)dy_n\cr &=\frac{1}{(n-1)!}\int_{\color{red}{y_n}=0}^1\left(\int_{[\color{red}{y_n},1]^{n-1}} \prod_{k=1}^{n-1}\frac{\ln(1-y_k)}{y_k}dy_1\ldots dy_{n-1}\right)\ln(1-y_n)dy_n \cr &=\frac{1}{(n-1)!}\int_{\color{red}{y_n}=0}^1\left(\int_{\color{red}{y_n}}^1 \frac{\ln(1-t)}{t}dt\right)^{n-1}\ln(1-y_n)dy_n \cr } $$ So, our first equivalent expression for $I_n$ is $$ I_n=\frac{1}{(n-1)!}\int_{0}^1\left(\int_{x}^1 \frac{\ln(1-t)}{t}dt\right)^{n-1}\ln(1-x)dx\tag 1 $$ Using integration by parts we have also $$ I_n=\frac{1}{n!}\int_{0}^1\left(\int_{x}^1 \frac{\ln(1-t)}{t}dt\right)^{n} dx\tag 2 $$ Noting that $\int_x^1\frac{\ln(1-t)}{t}dt={\rm Li}_2(x)-\frac{\pi^2}{6}$ we get $$\eqalign{ I_n &=\frac{1}{(n-1)!}\int_{0}^1\left({\rm Li}_2(x)-\frac{\pi^2}{6}\right)^{n-1}\ln(1-x) dx\cr &=\frac{1}{n!}\int_{0}^1\left({\rm Li}_2(x)-\frac{\pi^2}{6}\right)^{n} dx} \tag 3 $$ So, the question is reduced to evaluating the integral of powers of the dilogarithm. To my knowledge this is not known for powers larger than $2$.


Let $k$ be a positive integer, $s_1,\ldots,s_k$ positive integers with $s_1\geq 2$. The quantity $$ \tag{$\star$} \zeta(s_1,\ldots,s_k):=\sum_{n_1>\ldots>n_k\geq 1}\frac{1}{n_1^{s_1}\cdots n_k^{s_k}}\in\mathbb{R} $$ is called a multizeta value. The weight of $(\star)$ is $s_1+\ldots+s_k$ and the depth is $k$. Depth one multizeta values are just values of the Riemann zeta function. In weight up to seven, every multizeta value can be expressed in terms of Riemann zeta values, but this is conjecturally not true in weight eight and higher. For instance, it is conjectured that $\zeta(5,3)$ cannot be expressed in closed form in terms of the Riemann zeta function.

Every iterated integral of the form $$ \int_{\Delta}\frac{dt_1\ldots dt_n}{f_1(t_1)\ldots f_n(t_n)} $$ is a multizeta value, where $\Delta$ is the region $0\leq t_1\leq\ldots\leq t_n\leq 1$ and $f_i(t_i)$ is either $t_i$ or $1-t_i$, and every multizeta value can be expressed as such an integral. Moreover the conversion between integrals of this form and multizeta values is easy to obtain by expanding the integrand as a power series and integrating term by term.

In your integral $I_n$, substitute $y_i=x_1\ldots x_i$, and write each term $\log(1-y_i)$ as $-\int_0^{y_i} 1/(1-t_i)dt_i$. The resulting multiple integral can be divided up into regions depending on the relative ordings of the $y_i$ and $t_i$, and the integral over each region has the form $(\star)$. This means $I_n$ can be written as a linear combination of multizeta values.

This computation can be performed by computer (I'm using the Maple package HyperInt, written by Erik Panzer). After obtaining an expression for $I_n$, the package simplifies using known relations among multizeta values. For $n=4$ the result is $$ I_4=105-\frac{16}{7}\zeta(2)^3+6\zeta(3)^2-\frac{72}{5}\zeta(2)^2-30\zeta(3)-27\zeta(5)-\frac{65}{8}\zeta(7)+\frac{12}{5}\zeta(2)^2\zeta(3). $$ This has weight up to $7$, so can be expressed in terms of ordinary zeta values. The next integral $I_5$ has weight $9$ terms: $$ I_5=-945+{\frac {288}{7}}\,{\zeta(2)}^{3}-36\,{\zeta(3)}^{2}+108\, {\zeta(2)}^{2}+2\,{\zeta(3)}^{3}+210\,\zeta (3)+255\,\zeta (5)+{\frac {963}{8}}\,\zeta(7)+{\frac {8112}{875}}\,{\zeta(2)}^{4}-{\frac {24}{5}}\,\zeta (5,3)+{\frac {3299}{72}}\,\zeta (9)-36\,{\zeta (2)}^{2}\zeta (3)-{\frac {16}{7}}\,{\zeta (2)}^ {3}\zeta(3)-{\frac {66}{5}}\,{\zeta (2)}^{2}\zeta (5)-54\, \zeta (3)\zeta(5). $$ It's a huge mess, but notice the term $\zeta(5,3)$. So conjecturally $I_5$ cannot be written in terms of ordinary zeta values.

To summarize: there is a finite algorithm to express $I_n$ in terms of multizeta values, and for $n\geq 5$, we expect that $I_n$ cannot be expressed in terms of ordinary zeta values.