GOE. The most well known model of random matrices is probably the Gaussian Orthogonal Ensemble (GOE) of symmetric $n\times n$ random matrices, with Gaussian density proportional to $$X\mapsto\mathrm{e}^{-\tfrac{1}{2}\mathrm{Trace}(X^2)}.$$ The eigen-decomposition $X=ODO^\top$ and the change of variable $S\mapsto (O,D)$ give that the eigenvalues of GOE have density in $\mathbb{R}^n$ proportional to $$\mathrm{e}^{-\tfrac{1}{2}\sum_{i=1}^n\lambda_i^2}\prod_{i<j}|\lambda_i-\lambda_j|=\mathrm{e}^{-\Bigr(\tfrac{1}{2}\sum_{i=1}^n\lambda_i^2+\sum_{i<j}\log\frac{1}{|\lambda_i-\lambda_j|}\Bigr)}.$$ The term $\prod_{i<j}|\lambda_i-\lambda_j|$ is the absolute value of the determinant of the Jacobian of the change of variable. It expresses the way the Lebesgue measure on symmetric matrices is deformed by the change of variable. It vanishes when two eigenvalues are the same, revealing an eigenvalues repulsion. What is the simplest way to understand this repulsion phenomenon?

2×2 fact. It is actually a basic non-linear fact checkable on $2×2$ fixed matrices : $$X=\begin{pmatrix} x_1 & x_3\\x_3 & x_2\end{pmatrix},\quad x_1,x_2,x_3\in\mathbb{R}.$$ Namely, the eigenvalues $\lambda_1$ and $\lambda_2$, the roots of $\lambda^2-(x_1+x_2)\lambda+x_1x_2-x_3^2=0$, are $$\frac{x_1+x_2\pm\sqrt{(x_1+x_2)^2-4(x_1x_2-x_3^2)}}{2}$$ hence the eigenvalues spacing is given by $$|\lambda_1-\lambda_2|=\sqrt{(x_1-x_2)^2+4x_3^2}.$$ Now if $x_3=0$ then the matrix is diagonal and $|\lambda_1-\lambda_2|=|x_1-x_2|,$ while if $x_3\neq0$ then $|\lambda_1-\lambda_2|>2|x_3|>0$ which is precisely an eigenvalues repulsion phenomenon.

Surmise. Let us put now a Gaussian law on this $2×2$ matrix, making it GOE, with density $$\mathrm{e}^{-\tfrac{1}{2}\mathrm{Trace}\begin{pmatrix}x_1 & x_3\\x_3 & x_2\end{pmatrix}^2}=\mathrm{e}^{-\tfrac{1}{2}(x_1^2+x_2^2+2x_3^2)}.$$ Then the law of the spacing $S:=|\lambda_1-\lambda_2|$ can be computed as follows : $$S\overset{d}{=}\sqrt{(X_1-X_2)^2+4X_3^2},\quad\text{with } X_1,X_2,\sqrt{2}X_3\text{ iid }\mathcal{N}(0,1),$$ therefore
$$S^2\overset{d}{=}(X_1-X_2)^2+4X_3^2\sim\Gamma(\tfrac{1}{2},\tfrac{1}{4})*\Gamma(\tfrac{1}{2},\tfrac{1}{4})=\Gamma(1,\tfrac{1}{4})=\mathrm{Exp}(\tfrac{1}{4}),$$and thus the spacing $S$ has density proportional to $$s\mathrm{e}^{-\tfrac{1}{4}s^2}=\mathrm{e}^{-(\tfrac{1}{4}s^2-\log(s))}.$$ This density vanishes at the origin $s=0$, which corresponds to a logarithmic repulsion in the energy. This simple computation, made by Wigner, led him to his famous surmise (or conjecture) on the distribution of the eigenvalues spacings for high dimensional random matrices.

Non-linear effect and removable singularity. The squared spacing $S^2$ follows the exponential law of density proportional to $t\mapsto\mathrm{e^{-t/4}}$, which does not vanish at $t=0$, and this shows that the singularity at the origin of the density of the spacing $S$ can be removed by a non-linear change of variable. More generally, if for instance $t\mapsto f(t)$ is a density which may not vanish at $t=0$, then the change of variable $s=\sqrt{t}$ leads to a density proportional to $s\mapsto sf(s^2)$. In the case of a diagonal matrix, $x_3=0$, then $S^2=(X_1-X_2)^2\sim\Gamma(\tfrac{1}{2},\tfrac{1}{4})$, which has density proportional to $t\mapsto\mathrm{e}^{-t/4}/\sqrt{t}$, while $S$ has density proportional to $s\mapsto\mathrm{e}^{-s/4}$, which is not singular at $s=0$. What is important is the singularity or the non-singularity of the spacing $S$ because the squared spacing $S^2$ is not proportional to the coordinates of interest (the $\lambda$’s).

Jacobian. Let us check on the $2×2$ example that the absolute value of determinant of the Jacobian of the change of variable $X=ODO^\top\to (O,D)$ equals $|\lambda_1-\lambda_2|=\lambda_+-\lambda_-$ (here $n=2$). The matrix $X$ has $3=\frac{n^2-n}{2}+n$ degrees of freedom : $x_1,x_2,x_3$, while for $(O,D)$ we have two degrees of freedom $\lambda_\pm$, defining the diagonal matrix $D$, and another degree of freedom, say $\theta$, defining the $2×2$ orthogonal matrix $O$, which can be coded as $$O=\begin{pmatrix}\cos(\theta) & -\sin(\theta)\\\sin(\theta) & \cos(\theta)\end{pmatrix}.$$ Now $X=ODO^\top$ gives \begin{align*}x_1&=\lambda_+\cos(\theta)^2+\lambda_-\sin(\theta)^2\\x_2&=\lambda_+\sin(\theta)^2+\lambda_-\cos(\theta)^2\\x_3&=(\lambda_+-\lambda_-)\cos(\theta)\sin(\theta)\end{align*} and the determinant of the Jacobian is then \begin{align*}\det J&=\begin{vmatrix}\partial_{\lambda_+}x_1 & \partial_{\lambda_-}x_1 & \partial_{\theta}x_1\\\partial_{\lambda_+}x_2 & \partial_{\lambda_-}x_2 & \partial_{\theta}x_3\\\partial_{\lambda_+}x_3 & \partial_{\lambda_-}x_3 & \partial_{\theta}x_3\end{vmatrix}\\&=\begin{vmatrix}\cos(\theta)^2&\sin(\theta)^2&-2(\lambda_+-\lambda_-)\cos(\theta)\sin(\theta)\\\sin(\theta)^2&\cos(\theta)^2&2(\lambda_+-\lambda_-)\cos(\theta)\sin(\theta)\\\cos(\theta)\sin(\theta)&-\cos(\theta)\sin(\theta)&(\lambda_+-\lambda_-)(\cos(\theta)^2-\sin(\theta)^2)\end{vmatrix}\\&=\lambda_+-\lambda_-.\end{align*}

2×2 fact, non-Hermitian version. Let us consider now the fixed matrix $$X=\begin{pmatrix} x_1 & x_3\\x_4 & x_2\end{pmatrix},\quad x_1,x_2,x_3,x_4\in\mathbb{C}.$$ The eigenvalues $\lambda_1$ and $\lambda_2$, the roots in $\mathbb{C}$ of $\lambda^2-(x_1+x_2)\lambda+x_1x_2-x_3x_4=0$, are $$\frac{x_1+x_2\pm\sqrt{(x_1+x_2)^2-4(x_1x_2-x_3x_4)}}{2}$$ hence the eigenvalues spacing is given by $$|\lambda_1-\lambda_2|=\sqrt{|(x_1-x_2)^2+4x_3x_4|}.$$ Now if $x_3x_4=0$ then the matrix is tridiagonal and $|\lambda_1-\lambda_2|=|x_1-x_2|,$ while if $x_3x_4\neq0$ then the situation is more subtle than in the case of symmetric or Hermitian matrices.