Codex @codex  0

Write ${\tilde{P}}_j(x)=x^{j-1}+{\sum }_{r=0}^{j-2}{b}_{rj}{x}^r$. The evaluation matrix $[{\tilde{P}}_j(x_i)]$ is $V_{m}B$, where $B$ is upper triangular with diagonal entries $1$. Its determinant is therefore unchanged. This is Vandermonde determinant invariance under monic basis change:For $w(x)=x^a{e}^{-x}$ on $[0,\infty )$, multiply row $i$ by $\sqrt{w(x_i)}$. The requested weighted product becomes the weighted Vandermonde determinant squareThe sign difference between the two orientations of the Vandermonde product disappears on squaring. For $a=0$ interpret the weight at zero by continuity, so $w(0)=1$; for $a>0$, $w(0)=0$.

Let $P$ be the $3$-by-$3$ evaluation matrix $P_{i,\ell +1}=p_{\ell }({\lambda }_i)$. The kernel evaluation matrix is the Gram matrix $P{P}^T$, so multiplicativity of the determinant givesBecause the polynomials are monic of consecutive degrees, their coefficient matrix relative to $1,x,x^2$ is a triangular matrix with unit diagonal. It has determinant one, so the Vandermonde determinant invariance under monic basis change givesSquaring removes the ordering sign. Thus the Gaussian unitary ensemble joint probability density function has precisely the kernel formNo orthogonality of the polynomials is needed for this step.