Maximum of the Vandermonde determinant / minimum of the logarithmic energy

Write $V(a)$ for the determinant $\prod_{0\leq i<j\leq n-1} |a_i-a_j|$. Selberg's formula tells you that

$$\int_0^1 \cdots \int_0^1 V(a)^{2\beta} \prod_{i=0}^{n-1} da_i= n! \prod_{j=0}^{n-1} \frac{(\Gamma(1+j\beta))^2 \cdot \Gamma((j+1)\beta)} {\Gamma(2+(n+j-1)\beta)\cdot \Gamma(\beta)}=:A(n,\beta)$$

Thus the asymptotics you seek are given by $\lim_{\beta\to\infty} A(n,\beta)^{1/2\beta}$, which can be read from known asymptotics for the Gamma function. I did not try to perform the actual computation.

Remark: The constant $-2b^2$ is the maximum of the logarithmic energy $$\int \log |x-y| \mu(dx) \mu(dy) $$ over all probability measures supported on $[0,1]$. I am sure that this maximizer has been computed somewhere; Maybe it appears in Saff and Totik's book, which I do not have access to at the moment.

Since you know already that the optimal $a_i$ have $2a_i + 1 = x_i = \pm 1$ and the roots of $P'_{n-1}$, the calculation of $V_n$ comes down to the discriminant of $P'_{n-1}$, its leading coefficient, and its values at $\pm 1$, all of which are available in closed form via formulas for Jacobia polynomials (since $P'_{n-1}$ is a multiple of $P^{(1,1)}_{n-2}$). For the asymptotic growth of $\log M_n$, as ofer zeitouni suggests it is enough to find the maximum of the logarithmic energy $\int\!\!\int \log|x-y| \, d\mu(dx) \, d\mu(dy)$ over probability measures supported on $[0,1]$, and it is known that the optimal $\mu$ is the measure $\pi^{-1} dx/\sqrt{x-x^2}$ obtained from the uniform measure $d\theta$ on ${\bf R} / \pi{\bf Z}$ via $x = \cos^2 \theta$. The leading term is $\log M_n \sim -b n^2$ with $b=\log 2$, because the optimal measure on an interval of length $4$ has logarithmic energy zero so the average of the $n \choose 2$ terms for $(0,1)$ approaches $-\log 4$.

I am amazed that nobody noticed that Iosif asks for the calculation of the transfinite diameter of the interval $[0,1]$. This notion applies to arbitrary compact domains $K$ in ${\mathbb R}^n$ : $$d(K)=\lim_{k\rightarrow+\infty}\sup_{x_1,\ldots,x_k\in K}\left(\prod_{\alpha<\beta}|x_\alpha-x_\beta|\right)^{2/k(k-1)}.$$ In the particular case where $K\subset{\mathbb C}$ is simply connected, then the inverse of $d(K)$ is the conformal radius of ${\mathbb C}\setminus K$.