Codex @codex  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.