I would like to provide some details on the answers by **ofer zeitouni** and **Noam D. Elkies**.

By Selberg's formula, as presented in the answer by **ofer zeitouni**, one easily finds

$$(1)\qquad M_n=\lim_{\beta\to\infty} A(n,\beta)^{1/(2\beta)}=
\prod _{j=0}^{n-1} \frac{j^j (j+1)^{(j+1)/2}}{(j+n-1)^{(j+n-1)/2}}
$$
(with $0^0:=1$), which follows immediately from the observation that for any real $a$ and any real $c>0$
$$\Gamma(a+c\beta)^{1/(2\beta)}\sim(c\beta/e)^{c/2}
$$
as $\beta\to\infty$.

I asked for an upper bound on $M_n$, which would be asymptotic to $M_n$, thinking that an explicit expression for $M_n$ would not be possible. However, as is now clear from the answers by **ofer zeitouni** and **Noam D. Elkies**, such an expression is not so hard to obtain, and formula (1) presents this expression.

On the other hand, the relation between the logarithmic energy

$E_{n,a}:=2\sum_{0\le i<j\le n-1}\ln(a_j-a_i)$ and its ostensibly more general version $n^2E_\mu$ with $E_\mu:=\int\ln|x-y|\mu(dx) \mu(dy)$
for probability measures $\mu$ supported on $[0,1]$ seems unclear, I guess because of the singularities on the diagonal. Namely, one would expect that
$$(2)\qquad 2\ln M_n=\max_{0=a_0<\dots<a_{n-1}=1}E_{n,a}
\le n^2 \max_\mu E_\mu.
$$
However, it is not clear if $2E_{n,a}=2\sum_{0\le i<j\le n-1}\ln(a_j-a_i)$ could be written as
$n^2E_\mu=n^2\int\ln|x-y|\mu(dx) \mu(dy)$
for some probability measure $\mu$ on $[0,1]$.

In fact, quite surprisingly to me, the inequality in (2) is **false**, at least for large enough $n$. Indeed, it is not hard to show based on formula (1) that for some real $c\in(0,\infty)$
$$(3)\qquad M_n=(c+o(1))m_n,\quad\text{where}\quad m_n:=2^{-n^2} \sqrt{(n-1)!}\,(8e)^{n/2} n^{3/8}
$$
$$(4)\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad >>2^{-n^2}=\text{(?)}=\exp\Big\{\frac{n^2}2\,\max_\mu E_\mu\Big\}.
$$
The asymptotics $M_n=(c+o(1))m_n$ follows because
$$\ln\frac{M_{n+1}}{M_n}
=-n \ln2+\frac{n-1}{2}\,\ln (n-1)+\frac{n+1}{2}\,\ln (n+1)-\frac{2n-1}{2}\,\ln(2 n-1)$$
$$=d_n+O(1/n^2)$$
for large $n$, where
$d_n:=\frac{3}{8 n}-2n\ln2+\frac{1}{2}\ln n+\frac12(1+\ln2)$.

A curious corollary to (3)--(4) is that for each large enough natural $n$ there is some $a=(a_0,\dots,a_{n-1})$ with $0=a_0<\dots<a_{n-1}=1$ such that $2E_{n,a}$ cannot be approximated by (let alone written as) $n^2E_{\mu_k}$ for any sequence $(\mu_k)$ of probability measures on $[0,1]$.

Yet, it also follows from (3)--(4) that the logarithmic asymptotics $\ln M_n\sim\frac{n^2}2 \max_\mu E_\mu$ holds.

One can also see that $\ln(M_{n+1}/M_n)<d_n$ for $n\ge4$, and hence for each $k\ge4$ and all $n>k$ one has the upper bound $M_k\exp\sum_{j=k}^{n-1}d_j$ on $M_n$, which is asymptotic to $M_n$ as $n>k\to\infty$. In particular, for $n\ge4$ one has
$M_n<\tilde c_4 m_n$, where $m_n$ is as before and
$\tilde c_4:=\frac{512}{25} \sqrt{\frac{2}{15}} e^{(6 \gamma -43)/16}=0.631\dots$ and $\gamma$ is the Euler constant.