Spectral Density in the Gaussian Unitary Ensemble

In this post I perform numerical experiments on the spectral density in the Gaussian Unitary Ensemble (GUE).
This is the set of Hermitian \( N \times N \) matrices \( H \) with probability density \begin{equation*} P(H) dH = \mathcal{N} \prod_{ i < j} d \mathrm{Re} H_{ij} \prod_{ i < j} d \mathrm{Im} H_{ij}\ \prod_{ i} d H_{ii}\exp \left( -\frac{1}{2}\ \mathrm{tr}(H^2)\right) \end{equation*} where \( \mathcal{N} \) is the normalization constant. The spectral density \( \rho (\lambda) \) is defined by \begin{equation*} \rho(x) dx = \sum_{k=1}^N P(x \le \lambda_k \le x + dx) \end{equation*} with \( \lambda_k \) the \(N \) eigenvalues of the matrix \( H \). In chapter 6 of his book, Mehta calculates that \( \rho(x) = K_N(x,x) \) with \begin{equation*} K_N(x,y) = \frac{H_N(x) H_{N-1}(y) - H_N(y) H_{N-1}(x)}{x-y}e^{-\displaystyle\frac{x^2}{4} -\displaystyle\frac{y^2}{4}} \end{equation*} Thus \begin{equation}\label{eq:20151114c} \rho(x) = \left( H_{N-1}(x) H_N'(x) - H_{N}(x) H_{N-1}'(x) \right) e^{-\displaystyle\frac{x^2}{2}} \end{equation} Here \( H_N(x) \) are the Hermite polynomials. I now check formula \eqref{eq:20151114c} for \( N = 1, 2, 3 \) with a Monte Carlo simulation in Mathematica. I generate matrices from the GUE as follows: I draw matrices \( A \) and \( B \) whose elements \( A_{ij} \) and \( B_{ij} \) are independent normal random variables and take \begin{equation*} H = \frac{1}{2} \left( A + A^T \right) + \frac{i}{2} \left( B -B^T \right) \end{equation*} Then I calculate the \( N \) eigenvalues of \( H \). This is the code that I use for this.

Clear[randomsamples] (* This code calculates the eigenvalues of NSamples matrices from GUE. The size of the matrices is NN times NN *) 
randomsamples[NN_, NSamples_] := Table[A = RandomVariate[NormalDistribution[], {NN, NN}]; 
B = RandomVariate[NormalDistribution[], {NN, NN}]; 
H = 1/2 (A + Transpose[A]) + I / 2 (B - Transpose[B]); 
Eigenvalues[H], {k, NSamples}] // Flatten

N = 1
In this case \begin{equation*} \rho(x) = \frac{1}{\sqrt{2 \pi }} \exp \left(-\frac{x^2}{2}\right) \end{equation*} I produce a graph with the following commands  
NSamples = 10000; 
NN = 1; 
randomeigs = randomsamples[NN, NSamples]; 
fig1 = Histogram[randomeigs, {-3, 3, 0.1}, "PDF"]; 
fig2 = Plot[1/Sqrt[2 Pi] Exp[-x^2/2], {x, -3, 3}, PlotStyle -> Red]; 
Show[fig1, fig2]

N = 2
In this case \begin{equation*} \rho(x) = \frac{e^{-\frac{x^2}{2}} \left(x^2+1\right)}{2 \sqrt{2 \pi }} \end{equation*} I produce a graph with the following commands
NSamples = 50000; NN = 2; 
randomeigs = randomsamples[NN, NSamples]; 
fig1 = Histogram[randomeigs, {-4, 4, 0.1}, "PDF"]; 
fig2 = Plot[ (E^(-(x^2/2)) (1 + x^2))/(2 Sqrt[2 \[Pi]]), {x, -4, 4}, PlotStyle -> Red]; 
Show[fig1, fig2]

N = 3
In this case \begin{equation*} \rho(x) = \frac{e^{-\frac{x^2}{2}} \left(x^4+3\right)}{6 \sqrt{2 \pi }} \end{equation*} I produce a graph with the following commands
NSamples = 50000; NN = 3; 
randomeigs = randomsamples[NN, NSamples]; 
fig1 = Histogram[randomeigs, {-5, 5, 0.1}, "PDF"]; 
fig2 = Plot[ (E^(-(x^2/2)) (3 + x^4))/(6 Sqrt[2 \[Pi]]), {x, -5, 5}, PlotStyle -> Red]; Show[fig1, fig2]

The analytical expression \eqref{eq:20151114c} thus agrees well with the Monte Carlo simulation.
Warning: In this post, I have not paid attention to overall constants in expressions, I only made sure the normalizations were consistent in the Mathematica code.