Asymptotics of a Selberg-type integral Let $\Delta(s_1,s_2,\ldots,s_n) := \prod_{i<j}(s_i-s_j)^2$. Is there a standard way to estimate the decay of the Selberg-type integral
$$ I_n:=      \frac{1}{n!^2}\int_0^1 \int_0^1\cdots\int_0^1 \frac{\Delta(s_1,s_2,\ldots,s_n) \Delta(t_1,t_2,\ldots,t_n)}{\prod_{i,j}(1-s_i t_j)^2} d s_1\ldots d s_n d t_1 \ldots d t_n$$ as $n\to\infty$ ?
Checking numerically for $n\leq 25$, the decay seems to be like $e^{-2 n^2}$.
 A: You can follow large deviation-type estimates from random matrix theory, starting with Ben Arous & Guionnet's paper: https://link.springer.com/article/10.1007/s004400050119
You will eventually obtain that your Selberg-type integral behaves like 
$$
\frac{1}{(n!)^2}e^{-n^2 E_*}
$$
where $E_*$ is the minimum of the functional 
\begin{multline}
E(\mu,\nu)=\iint \log\frac1{|x-y|}d\mu(x)d\mu(y)\\+\iint \log\frac1{|x-y|}d\nu(x)d\nu(y)\\-2\iint\log\frac1{|1-xy|}d\mu(x)d\nu(y)
\end{multline}
where $(\mu,\nu)$ ranges over pairs of Borel probability measures supported in $[0,1]$. I guess you can work out the existence and unicity of a minimizer for $E$ and obtain Euler-Lagrange equations to characterize this minimiser and, hopefully, compute explicitly the minimal value $E_*$.
There is quite some technical work to do to fill the gaps, and sorry for the self-advertisement, but if you follow this approach let me tell you that I had to struggle with a similar two type particle problem in the paper http://www.worldscientific.com/doi/abs/10.1142/S2010326312500165 written with Arno Kuijlaars. I hope this can help. 
By the way, minimizing $E$ is referred as to a vector equilibrium problem in potential theory. 
ADDENDUM: After thinking further, it seems more natural to make the changes of variables $s_i=1/x_i$, which leads to the same asymptotics but with $E_*$ the minimum of 
\begin{multline}
\tilde E(\mu,\nu)=\iint \log\frac1{|x-y|}d(\mu-\nu)(x)d(\mu-\nu)(y)+\int_1^\infty \log(x)d\mu(x)
\end{multline}
where $\mu,\nu$ are Borel probability measures living on $[1,\infty)$ and $[0,1]$ respectively. 
A: This is very late to answer this question, but let me mention a rather different approach.
You can use the Schur function expansion
$$\frac{1}{\prod_{ij}(1-s_it_j)^2}=\sum_{\lambda,\mu}S_\lambda(s)S_\lambda(t)S_\mu(s)S_\mu(t),$$
where the infinite sums are over partitions, to get
$$I_n=\frac{1}{n!^2} \sum_{\lambda,\mu}(g_{\lambda,\mu})^2,$$ with
$$g_{\lambda,\mu}=\int_0^1 \Delta(s) S_\lambda(s)S_\mu(s)ds.$$
The above integral is known (see L.K. Hua, Harmonic analysis of functions of several complex variables in the complex domains, Translations of mathematical monographs, Vol. 6) to be
$$g_{\lambda,\mu}=n!\frac{\prod_{1\le i<j\le n}(\lambda_i-\lambda_j+j-i)(\mu_i-\mu_j+j-i)}{\prod_{1\le k,r\le n}(2n+\lambda_k+\mu_r-k-r+1)}.$$
I am not sure how to extract the asymptotics from the above expression, though.
