Demo Notebook illustrating  Type I class matrices.

"A Class of Parameter Dependent Comuting Matrices"

Sriram Shastry 
Physics Department, UCSC, Santa Cruz, CA 95064
sriram@physics.ucsc.edu
Cond- mat  0501502 January 20, 2005

[Version 1.0 of 20 Jan 2005]

This notebook contains a numerical example in d=30 where all the assignable variables of the matrix are given random values. The constraints are used to compute the remaining matrix elements. 
We construct the S matrix first, then the α matrix, followed by the β matrix. In this demo program, you can see that the α and β actually commute, the spaghetti of eigenvalues of α as well as that of β shows no signs of level repulsion. Try using d= 6 or 10 to get a closer look at some nice level crossings. At the end of the notebook, one sees that the energy level distribution histograms continue to have peaks near the origin, as one expects from Poisson statistics. 

Feel free to  play with the various parameters, e.g. the size of the random numbers, to generate  other interesting  cases!

d :=  30;

R[i_, j_] :=  - R[j, i] /; i > j
eq[i_, j_, k_, l_] := 
      R[i, j] R[k, l] - R[i, k] R[j, l] + R[i, l] R[j, k]
varset := Flatten[Table[R[i, j], {i, 3, d}, {j, i + 1, d}]]
expressionset := 
        Flatten[Table[
            eq[1, 2, i, j], {i, 3, d}, {j, i + 1, 
              d}]];
\(\(ll = \((d - 2)\) \((d - 3)\)/2;\)

\)

SeedRandom[1950]
rr := Random[Real, {-2, 2}, 7]
rroffdiagonal := 
      Random[Real, {-1/Sqrt[d], 1/Sqrt[d]}, 7]

Random[Real,{-1,1},7] generates a random number in the interval -1 to 1 with 7 digits. The idea of rr and rroffdiagonal is that these are random numbers on scale 1 and 1/Sqrt[d] respectively, this scaling assures that the α matrix will have eigenvalues of scale unity.

givenvars = 
      Block[{}, Do[R[1, i] = 1/rroffdiagonal, {i, 2, d}]; 
        Do[R[2, j] = 1/rroffdiagonal, {j, 3, d}]]
solver = 
        Solve[Table[expressionset[[j]], {j, 1, ll}] == 
            Table[0, {j, 1, ll}], varset];
Do[S[i, j] = (1/R[i, j]) /. solver[[1]], {i, 1, d - 1}, {j, 
        i + 1, d}]

Thus the S matrix has now  been fixed...next we move on to  fixing the μ type equations.

S[i_, i_] = 0; S[i_, j_] := -S[j, i] /; j < i
p[i_, j_] := (A[i] - A[j])/(a[i] - a[j])
Y[i_, j_] := 
      p[i, j] - 
        1/(S[i, j] ( a[i] - a[j])) Sum[
            S[i, l] S[l, j] ( a[l] - a[j]) , {l, 1, d}]
xi[i_, j_] := 
      p[i, j] (b[i] - b[j]) + 
        1/(S[i, j] (a[i] - a[j])) 
          Sum[ S[i, l] S[l, 
                j] (b[i] a[j] - b[j] a[i]  + b[j] a[l] - b[l] a[j] + 
                  b[l] a[i] - b[i] a[l]), {l, 1, d}]
mueq[j_] := Y[1, 2] - Y[1, j] - S[1, 2] S[1, j]/S[2, j]
muvarset := Table[A[j], {j, 3, d}]
solvemu = Solve[
            Table[mueq[j], {j, 3, d}] == Table[0, {j, 3, d}], 
            muvarset][[1]];

specifya = {Join[
              Table[a[i] -> rr, {i, 1, d}], {A[1] -> rr, 
                A[2] -> rr}]}[[1]];
α[x_] =  
        DiagonalMatrix[Table[a[i]  +  x A[i] /. solvemu, {i, 1, d}]] +  
            x Table[(a[i] - a[j]) S[i, j], {i, 1, d}, {j, 1, d}] /. 
          specifya;

The above commands gives a  realizations of the alpha matrix where the diagonal entries a[j] are chosen randomly, and
so also  A[1],A[2].

dontshow = {DisplayFunction -> Identity};
doshow = {DisplayFunction -> $DisplayFunction};
getmultiple := << Graphics`Master`
pj = {PlotJoined -> True};
mlp[list_] := 
      Block[{}, 
        Do[gr[i] = ListPlot[list[[i]], dontshow, pj], {i, 1, 
            Length[list]}]; 
Show[
          Table[gr[i], {i, 1, Length[list]}], doshow]]
evgrid[a_, xmin_, xmax_, nx_] := 
      Block[{y, dx, res1}, dx = (xmax - xmin)/nx; 
        y[j_] = xmin + (j - 1) dx; 
        res1 = d Table[
              Sort[Eigenvalues[a[y[r]]]], {r, 1, 
                nx}]; 
lers = Length[res1]; 
        Do[g[n] = Table[{y[j], res1[[j, n]]}, {j, 1, nx}], {n, 1, d}]; 
        Return[Table[g[n], {n, 1, d}]]]

"mlp" is my personal  command for multiple list plotting.

mlp[evgrid[α, - .1,  .1, 1000]]

Now we get about β(x) that commutes with α(x). We have the d-1 equations xi(1,j)=0 that can fix B[j] for j=2,..d provided we specify b[j] and B[1]

specifyb = {Join[
              Table[b[i] -> rr, {i, 1, d}], {B[1] -> 
                  rr}]}[[1]];
Do[
        B[j] = B[1] +  xi[j, 1] //. Join[specifya, specifyb, solvemu], {j, 
          2, d}];
β[x_] = 
        DiagonalMatrix[
            Table[( b[i] /. specifyb ) + ( 
                  x B[i] /.  specifyb), {i, 1, d}]] +  
          x Table[
              S[i, j] (b[i] - b[j]) /. specifyb, {i, 1, d}, {j, 1, 
                d}];

Let us check if α and β commute.

Chop[
      Expand[α[x] . β[x] - β[x] . α[x]]]

This shows that these matrices actually do commute for all values of x.

mlp[evgrid[β, -1, 1, 1000]]

Notice the profusion of level crossings. Take a smaller value of d or range to see more details.

<< Statistics`DataManipulation`
<< Graphics`Graphics`
levelspectrum[a_] := 
      Block[{}, t0 = Sort[Eigenvalues[a]]; t1 = Min[t0]; t2 = Max[t0]; 
        dx = (t2 - t1)/Sqrt[ d]; t3 = BinCounts[t0, {t1, t2, dx}]; 
        midpoints = Table[x + dx/2, {x, t1, t2, dx}]; 
        Return[BarChart[
            Table[{t3[[i]], midpoints[[i]]}, {i, 1, 
                Length[t3]}]]]]

levelspacer[a_] := 
      Block[{}, t1 = Sort[Eigenvalues[a]]; scale = (Max[t1] - Min[t1])/d; 
        t2 = 1/scale Table[
              t1[[i]] - t1[[i - 1]], {i, 2, Length[t1]}]; 
        Return[Sort[t2]]]
levelstatistics[a_, dx_] := 
      Block[{}, temp1 = levelspacer[a]; temp2 = Max[temp1]; 
        temp3 = BinCounts[temp1, {0, temp2, dx}]; 
        midpoints = Table[x + dx/2, {x, 0, temp2, dx}]; 
        Return[BarChart[
            Table[{temp3[[i]], midpoints[[i]]}, {i, 1, 
                Length[temp3]}]]]]

levelstatistics[α[1], 1]

levelstatistics[β[1], 1]

levelspectrum[α[1. ]]