$\newcommand{\trace}{\operatorname{trace}}$
The result below mentions a reasonably improved inequality.
Let $m = \frac{\trace(A)}{n}$, and $s^2= \frac{\trace(A^2)}{n}-m^2$. Then, Wolkowicz and Styan (Linear Algebra and its Applications, 29:471-508, 1980), show that
\begin{equation*}
\lambda_1 \ge \frac{\det(A)}{(m+s/\sqrt{n-1})^{n-1}}
\end{equation*}
Remark: As per the notation in the OP, $\lambda_1$ is the smallest eigenvalue---usually the literature uses $\lambda_1$ to be largest.
Thus, we obtain the upper bound
\begin{equation*}
\lambda_2\lambda_3\cdots\lambda_n \le \left(m+ \frac{s}{\sqrt{n-1}}\right)^{n-1}.
\end{equation*}
This bound is tight. Consider for example,
If $A= \text{Diag}(1,2,2,\ldots,2)$, then the lhs is $2^{n-1}$, $m=2-1/n$, and $s^2
= 1/n-1/n^2$, so that $s/\sqrt{n-1} = 1/n$. Thus, the bound on the rhs is tight.
With $M := \max_{i,j}|a_{ij}|$, we see that $m \le M$ and $s^2 \le nM^2 - m^2$, which leads to an upper bound in terms of $M$ as desired
\begin{equation*}
\lambda_2\lambda_3\cdots\lambda_n \le \left(M+ \sqrt{\frac{nM^2-m^2}{n-1}}\right)^{n-1} < \left(M + M\sqrt{\frac{n}{n-1}} \right)^{n-1},
\end{equation*}
which is better than the bound mentioned in the post (though we lost a bit by deleting $m^2$).
The limit is positive but well below $1$; I show that the limit is in $(0.09448, 0.2646)$, and outline how to approximate it much more closely than this.
The average of $\mu$, call it ${\rm E}(\mu)$,
is twice the average number of maximal all-0 polyominos;
so the limit of ${\rm E}(\mu) / n^2$ is twice the expected number of
maximal all-0 polyominos per unit area of a large square.
That limit is the special case $p=q=1/2$ of
the expected number, call it $\epsilon(p)$,
of maximal all-0 polyominos per unit area
when each entry is independently 0 or 1 with probability $p,q$ respectively
($p+q = 1$). We can approximate $\epsilon(p)$ by writing it as a sum
$\sum_P \epsilon_P(p)$ over all possible polyomino shapes $P$,
and calculating partial sums over all $P$ of size at most $k$.
The contribution $\epsilon_P(p)$ is $p^{|P|} q^{|\delta(P)|}$,
where $|P|$ is the size of $P$, and $|\delta(P)|$ is the number of
cells at distance $1$ from $P$. (We do not identify different orientations;
for example, for $k=4$ there is one square $P$, two straight ones,
four S/Z shaped tetrominos, four T-shaped, and eight that are L-shaped.)
So we seek $\sum_P p^{|P|} q^{|\delta(P)|}$. Note that
$\sum_P |P| p^{|P|} q^{|\delta(P)|} = p$ because this is just
the expected number of 0 entries per unit area.
The analysis so far works in any dimension; for example, in dimension $1$
there is a unique $P$ of each size $k \geq 1$, and $|\delta(P)| = 2$ always,
so $\epsilon(p) = \sum_{k=1}^\infty p^k q^2 = p q^2/(1-p) = p q$
(and indeed $\sum_{k=1}^\infty k p^k q^2 = p q^2/(1-p)^2 = p$).
For example, $\epsilon(1/2) = 1/4$, and the expected number of
maximal all-0 or all-1 strings in a random bitstring of length $n$
is asymptotically $2\epsilon(1/2) n = n/2$.
(Exercise [noted in my comment]: in fact that expected number is
exactly $(n+1)/2$. What happens for arbitrary $p$?)
In dimension $2$, the sum of $\epsilon_P(p)$ over $|P| \leq 4$ is
$$
p q^4 + 2 p^2 q^6 + p^3 (2 q^8 + 4 q^7)
+ p^4 (q^8 + 2 q^{10} + 4 q^8 + 4 q^8 + 8 q^9),
$$
and the sum of $|P| \epsilon_p(P)$ up to $|P| \leq 4$
is obtained from by multiplying the $k$-th term by $k$ ($k=1,2,3,4$),
and comes to $p - 315 p^5 + 1824 p^6 - 4848 p^7 + \cdots$.
(The coefficient $315$ seems to come from the number of rooted 5-ominos,
see https://oeis.org/A048664 .)
For $p=1/2$ these sums are only $387/2^{13}$ and $612/2^{13}$
respectively. It follows that the $|P| \geq 5$ terms must contribute
$1/2 - 612/2^{13}$ to the sum of $|P| \epsilon_p(P)$,
and thus can contribute at most $1/5$ of that to the sum of $\epsilon_p(P)$.
We conclude that
$$
\frac{387}{2^{13}} < \epsilon(1/2) <
\frac{387}{2^{13}} + \frac15 \left( \frac12 - \frac{612}{2^{13}} \right)
= \frac{5419}{40960};
$$
numerically these bounds are $0.04724+$ and just under $0.1323$.
Hence the answer to the OP's question is between $0.09448$ and $0.2646$,
as claimed.
Polyominos have been tabulated well beyond $k=4$.
It must be feasible to compute $\sum_{|P| \leq k} \epsilon_P(1/2)$ and
$\sum_{|P| \leq k} |P| \epsilon_P(1/2)$ for $k$ large enough to obtain
a reasonably close approximation to $\epsilon(1/2)$, and thus to
the value $2\epsilon(1/2)$ of the limit that the OP seeks.
Best Answer
As noted in Will's comment above, it's easy to compute the expected square of the determinant. More precisely, we have $$E(\det M^2)=n! \prod \frac{k_i (k_i+1)}{3}.$$
Let $M'$ be formed from $M$ by dividing each row by $\left(\frac{k_i(k_i+1)}{3}\right)^{1/2}$. Now each entry has mean $0$ and variance $1$, and furthermore the entries are bounded. The determinants of such matrices as the size of the matrix tends to infinity have been well studied, by Girko, Tao and Vu, and Nguyen and Vu, among others.
For example, it follows from Theorem 1.1 in the Nguyen and Vu paper linked above that $\log |det(M')|$ is asymptotically normal with mean $\frac{1}{2} \log((n-1)!)$ and variance $\log n$. Taking this and rescaling back to $M$, we have that with probability tending to $1$ as $n$ tends to infinity that $$ det M^2 = \ n^{-1+o(1)} E(det M^2).$$
In particular, the squared determinant is almost surely concentrated in a short interval which does not contain its expectation!