How to compute $I(n)=\int_{0}^1x_1^a\int_{x_1}^1x_2^a\cdots\int_{x_{n-1}}^{1}x_n^adx_n\cdots dx_2dx_1$? Let $a$ be a positive integer. For $n\ge 1$, how to compute
$$I(n):=\int_{0}^1x_1^a\int_{x_1}^1x_2^a\cdots\int_{x_{n-1}}^{1}x_n^adx_n\cdots dx_2dx_1?$$

I tried first two terms, and found $I(1)=\frac{1}{a+1}$, $I(2)=\frac{1}{2(a+1)^2}$. But I don't see how to compute $I(n)$ by $I(n-1)$ inductively.
 A: Method 1. We have
$$ I(n) = \int_{0 \leq x_1 \leq \cdots \leq x_n \leq 1} x_1^a \cdots x_n^a \, \mathrm{d}x_1\cdots\mathrm{d}x_n. $$
Now by symmetry, this is simply
$$ I(n) = \frac{1}{n!} \left( \int_{0}^{1} x^a \, \mathrm{d}x \right)^n = \frac{1}{n!(a+1)^n}. $$


*

*Addendum. More specifically, for each permutation $\sigma$ on $[n]=\{1,\cdots,n\}$, define
$$\Delta(\sigma) = \{ (x_1, \cdots, x_n) \in \mathbb{R}^n : 0 \leq x_{\sigma(1)} \leq \cdots \leq x_{\sigma(n)} \leq 1\}$$
If $S_n$ denotes the set of all permutations on $[n]$, the sets $\Delta(\sigma)$ are non-overlapping for different $\sigma$'s in $S_n$ and $\bigcup_{\sigma \in S_n} \Delta(\sigma) = [0, 1]^n$. Moreover, by symmetry, we may as well write
$$ \forall \sigma \in S_n, \quad I(n) = \int_{\Delta(\sigma)} x_1^a \cdots x_n^a \, \mathrm{d}x_1\cdots\mathrm{d}x_n. $$
(This follows from the substitution $(x_1, \cdots, x_n) \mapsto (x_{\sigma(1)}, \cdots, x_{\sigma(n)})$.) Therefore
$$ n!I(n)
= \sum_{\sigma \in S_n} I(n)
= \int_{\bigcup \Delta(\sigma)} x_1^a \cdots x_n^a \, \mathrm{d}x_1\cdots\mathrm{d}x_n
= \int_{[0,1]^n} x_1^a \cdots x_n^a \, \mathrm{d}x_1\cdots\mathrm{d}x_n. $$
Now the latter integral factors out into product of $\int_{0}^{1} x_i^a \, \mathrm{d}x_i$'s, proving the desired claim.

Method 2. Alternatively, define
$$ y(t) = 1 + \sum_{n=1}^{\infty} I(n) t^{n(a+1)} = 1 + \sum_{n=1}^{\infty} \int_{0<x_1<\cdots<x_n<t} x_1^a \cdots x_n^a \, \mathrm{d}x_1 \cdots \mathrm{d}x_n. $$
This solves the following integral equation
$$ y(t) = 1 + \int_{0}^{t} s^a y(s) \, \mathrm{d}s. $$
This is equivalent to $y(0) = 0$ and $y'(t) = t^a y(t)$, and so, $y(t) = e^{t^{a+1}/(a+1)}$. From this, we easily read out the value of $I(n)$.
