Closed-Form solution for nested integrals of this polynomial? I was wondering whether there is a closed-form solution for this (nested) integral:
$$
\int_{-1}^{1}\int_{t_{0}}^{1}\int_{t_{1}}^{1}...\int_{t_{a-2}}^{1}\prod_{\begin{array}{c}
i<j\\
j=\{0,..,a-1\}\\
i=\{0,..,a-1\}
\end{array}}\left(t_{i}-t_{j}\right)^{4}dt_{a-1}dt_{a-2}...dt_{0}
$$
These are the results I get for $a=2$ and $a=3$:
$a=2$:
$$
\int _{-1}^1\int _{t_0}^1(t_0-t_1)^4 dt_1 dt_0 = \frac{32}{15}
$$
$a=3$:
$$
\int _{-1}^1\int _{t_0}^1\int _{t_1}^1(t_0-t_1)^4 (t_0-t_2)^4 (t_1-t_2)^4dt_2dt_1dt_0 = \frac{8192}{33075}
$$
Is there a known closed-form solution $\forall a$?
Edit: As noted by  @Steven Stadnicki in the comments, the function inside the integrals can also be written as
$$
\prod_{\begin{array}{c}
i\ne j\\
j=\{0,..,a-1\}\\
i=\{0,..,a-1\}
\end{array}}(t_{i}-t_{j})^{2}
$$
Thanks!
 A: I'm using (verbatim) the easy-to-read paper [1] to show that
$$J_n:=\idotsint\limits_{-1<t_1<\ldots<t_n<1}\prod_{1\leqslant j<i\leqslant n}(t_i-t_j)^4\,dt_1\ldots dt_n\color{blue}{=2^{n(2n-1)}\prod_{0<j<i<2n}\frac{i-j}{i+j}}.\tag{*}\label{result}$$
[1] N.G. de Bruijn, On some multiple integrals involving determinants, 1955
Observe that $J_n=\frac{1}{n!}\idotsint_{[-1,1]^n}=\frac{2^{n(2n-1)}}{n!}\idotsint_{[0,1]^n}$ (the first equality follows from the symmetry of the integrand; the second one is obtained after substituting $t_i=2x_i-1$ and renaming $x_i$ back to $t_i$).
Now the relevant result from the paper (see section $7$) is as follows. For $1\leqslant i,j\leqslant2n$, define $$F_{i,j}(t_1,\ldots,t_n)=\begin{cases}\varphi_i(t_k),&j=2k-1\\\psi_i(t_k),&j=2k\end{cases},\quad G_{i,j}=\int_\Omega\varphi_i(t)\psi_j(t)\,dt,$$ where $\varphi_k,\psi_k : \Omega\to\mathbb{R}$ are good enough (for all the integrals to exist), and let $F(t_1,\ldots,t_n)$ and $G$ be the corresponding $2n\times2n$ matrices. Then $$\idotsint_{\Omega^n}\det F(t_1,\ldots,t_n)\,dt_1\ldots dt_n=2^n{n!}\operatorname{Pf}\widehat{G},$$ where $\operatorname{Pf}\widehat{G}$ is the Pfaffian of $\widehat{G}=\frac12(G-G^\mathsf{T})$.
In our case, we choose $\Omega=[0,1]$, $\varphi_k(t)=t^{k-1}$ and $\psi_k(t)=\varphi_k'(t)$; from the answer to a recent question of yours (see also this one by myself), we know that $\det F(t_1,\ldots,t_n)$ is precisely our integrand in $J_n$. Thus, $J_n=2^{2n^2}\operatorname{Pf}\widehat{G}$ where $G_{1,1}=\widehat{G}_{1,1}=0$ and, otherwise, $$G_{i,j}=(j-1)\int_0^1 t^{i+j-3}\,dt=\frac{j-1}{i+j-2}\implies 2\widehat{G}_{i,j}=\frac{j-i}{i+j-2}.$$
So, it remains to compute $\operatorname{Pf}\widehat{G}$ or its square $\det\widehat{G}$. Now the relevant result from the paper is $$\det_{1\leqslant i,j\leqslant 2n}\frac{x_i-x_j}{x_i+x_j}=\prod_{1\leqslant i<j\leqslant 2n}\left(\frac{x_i-x_j}{x_i+x_j}\right)^2$$ (stated near the beginning of section $9$; see this question). This leads to \eqref{result} directly.
