Solution (source code)

= Solution

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> $PP^T$, so <multiplicativity of the determinant> gives
$$
\det[K_3(\lambda_i,\lambda_j)]_{i,j=1}^3
=\det(PP^T)=(\det P)^2.
$$
Because 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> gives
$$
\det P=\det\begin{pmatrix}
1&\lambda_1&\lambda_1^2\\
1&\lambda_2&\lambda_2^2\\
1&\lambda_3&\lambda_3^2
\end{pmatrix}
=\prod_{i<j}(\lambda_j-\lambda_i).
$$
Squaring removes the ordering sign. Thus the <Gaussian unitary ensemble> <joint probability density function> has precisely the kernel form
$$
\boxed{\rho_3(\boldsymbol\lambda)
=c_3\det[K_3(\lambda_i,\lambda_j)]_{i,j=1}^3
\exp\!\left(-\frac12\sum_i\lambda_i^2\right).}
$$
No <orthogonality> of the <polynomials> is needed for this step.