How to compute the following integral in $n$ variables?
Solution 1:
In this answer, I will provide a method that uses iterated partial summation on the alternating Harmonic numbers that will allow you to explicitly evaluate either of these integrals for any $n$.
Start by expanding the series $$\frac{1}{1-\prod_{k}x_{k}}=1+\prod_{k}x_{k}+\prod_{k}x_{k}^{2}+\cdots$$ so that $$I_{n}=\sum_{j=0}^{\infty}\left(\int_{0}^{1}\frac{1-x}{1+x}x^{j}dx\right)^{n}=\sum_{j=0}^\infty r_j^n.$$ Next we will rewrite $r_j$ in a more manageable form. Expanding the power series we have that $$r_j=\int_{0}^{1}\frac{1-x}{1+x}x^{j}dx=\sum_{k=0}^{\infty}(-1)^{k}\int_{0}^{1}(1-x)x^{j+k}dx,$$ and by the Beta function identity $$\text{B}(x,y)=\frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}$$ we obtain $$r_j=\sum_{k=0}^{\infty}(-1)^{k}\frac{1}{(k+j+1)(k+j+2)}.$$ The terms in this series may be split by partial fractions yielding $$(-1)^{j}\sum_{k=0}^{\infty}\frac{(-1)^{k+j}}{k+j+1}+(-1)^{j}\sum_{k=0}^{\infty}\frac{(-1)^{k+j+1}}{k+j+2}$$ and since $\sum_{k=0}^\infty (-1)^k/(k+1)=\log2$ we have $$r_j=(-1)^{j}\left(2\log2-H_j'-H_{j+1}'\right)$$ where $$H_j'=\sum_{k=1}^{j}\frac{(-1)^{k-1}}{k}$$ is the $j^{th}$ alternating harmonic number. To evaluate the series, we need only compute the partial sums of $$\sum_{j=0}^m (-1)^{nj}\left(2\log2-H_j'-H_{j+1}'\right)^n$$ for any $n,m$. This can be done explicitly using the technique of iterated partial summation as shown in this answer. This will yield an exact answer for any $n$.
Solution 2:
Since all the $a_k \lt 1$ (except at one point), I would try expanding $$\frac{1}{1-\prod_k^n a_k}=1+\left(\prod_k^n a_k \right)+\left(\prod_k^n a_k \right)^2+\cdots$$ (Geometric series) This will allow you to turn your problem into a sum of integrals in which, in every term of the sum, the $a_k$'s appear multiplicatively and thus can be integrated independently via Fubini:
$$\sum_{i=0}^\infty \int_0^1 \cdots \int_0^1 \prod_{k=1}^n \left( \frac{1-a_k}{1+a_k} (a_k)^i\right) $$. Which is equal to
$$\sum_{i=0}^\infty \left( \int_0^1 \left( \frac{1-a_k}{1+a_k} (a_k)^i\right) \right)^n $$
I can't assure this will work, but it seems to me that the problem has been reduced to computing one integral (which doesn't look impassable) and a series.
Hope this helps and let me know if it does.