Random Matrix Theory Hints
(developed with Piet Brouwer)

We’re generating M members of the GOE ensemble of NxN matrices, and calculating the difference between the eigenvalues in the middle of the range. Let’s demonstrate how using a small sample.

M=200;
Size=4;

To generate random symmetric matrices, we generate a matrix with entries drawn from a normal distribution, and then symmetrize it.

Mat = RandomVariate[NormalDistribution[],{Size,Size}]
MatSym = Mat + Transpose[Mat]
MatrixForm[MatSym]

Generate M random SizexSize arrays of Gaussian (normal) random numbers;symmetrize them; find their sorted eigenvalues,keep a list of the difference between the two middle ones. You can do this with a For loop (slow and inelegant)

diffs=Table[0,{i,1,M}];
M11=Table[0,{i,1,M}];
For[i=1,i≤M,i++,{Mat=RandomVariate[NormalDistribution[],{Size,Size}];
MatSym=Mat+Transpose[Mat];
eigs=Sort[Eigenvalues[MatSym]];
diffs[[i]]=eigs[[Size/2+1]]-eigs[[Size/2]]};
  M11[[i]]=MatSym[[1]][[1]];]

You can also do this with functional programming: it' s considered more elegant in Mathematica world

Symmetrized[m_]:=m+Transpose[m];
EigSorted[m_]:=Sort[Eigenvalues[m]]
EigDiffs[eigs_]:=eigs[[Size/2+1]]-eigs[[Size/2]]
diffs=Map[EigDiffs[EigSorted[Symmetrized[#]]]&,RandomVariate[NormalDistribution[],{M,Size,Size}]];
M11 = Map[Symmetrized[#][[1,1]]&,Table[Random[NormalDistribution[]],{m,1,M},{i,1,Size},{j,1,Size}]];

Or, perhaps, with tables :

Mats = RandomVariate[NormalDistribution[],{M,Size,Size}];
MatsSym = Table[Mats[[i]]+Transpose[Mats[[i]]],{i,1,M}];
M11 = Table[MatsSym[[i]][[1,1]],{i,1,M}];
eigs =Table[Sort[Eigenvalues[MatsSym[[i]]]],{i,1,M}];
diffs = Table[eigs[[i]][[Size/2+1]]-eigs[[i]][[Size/2]],{i,1,M}];

Divide out by mean value of the splittings

diffAve=Mean[diffs]

Histogram generates a histogram of diffs. ‘Automatic’ tells it to choose the number of bins. “Probability” normalizes the integral to one.

Histogram[diffs/diffAve,Automatic,"PDF"]

Let’s demonstrate how to compare simulation with theory by using a histogram for the diagonal element M11. First plot the histogram for M11; store as "histGauss" for later plotting with theory

histGauss=Histogram[M11,Automatic, "PDF"]

(*Now,plot expected Gaussian fit*)
sigma=2;
theoryGauss=Plot[(1/Sqrt[2 Pi sigma^2]) Exp[-x^2/(2 sigma^2)],{x,-5,5}]

(*and show the comparison*)

Show[histGauss,theoryGauss]

We generate random +-1 matrices using RandomInteger (which produces 0/1, which we double and subtract one). Note that Mathematica subtracts one ‘componentwise’. For example...

MatPM=RandomInteger[{0,1},{Size,Size}]*2-1

This is an inelegant way of symmetrizing ...

MatPMSym=Table[If[i≥j,MatPM[[i]][[j]],MatPM[[j]][[i]]],{i,1,Size},{j,1,Size}];

MatrixForm[MatPMSym]