Spin ring:  N=6, spin=1/2

Heisenberg model,
complete diagonalization

©  Apl. Prof. Dr. Jürgen Schnack, University of Osnbrueck, Germany

Initialization

Set -2 J = 1; (AF)
Basis is one position longer to simulate NS+1=1.

NS = 6;
spin = (1/2);
\(\(\(DimH = \((2*spin + 1)\)\^NS;\)\(
\)
    \)\)
Basis = Table[{0, 0, 0, 0, 0, 0, 0}, {n, 1, DimH}];

Create basis

For[i = 2, 
        i < DimH + 1, i++, 
For[j = 1, 
          j < NS + 1, j++, 
Basis[[i, j]] = 
            Basis[[i - 1, j]]]; 
Basis[[i, 1]] = 
          Basis[[i, 1]] + 1; 
j = 
          1; 
While[
          j < NS, 
If[
            Basis[[i, j]] == (2*spin + 1), 
            Basis[[i, j]] = 0; 
            Basis[[i, j + 1]] = Basis[[i, j + 1]] + 1; j++, 
            j = NS]]
];
For[i = 1, 
        i < DimH + 1, i++, 
Basis[[i, NS + 1]] = 
          Basis[[i, 1]]];

MatrixForm[Basis]

Create full Hamilton matrix

Diagonal = Ising model

HMatrix = Table[0.0, {i, 1, DimH}, {j, 1, DimH}];
\(\(For[i = 1, 
        i < DimH + 1, \(i++\), 
For[j = 1, 
          j < DimH + 
              1, \(j++\), 
		HMatrix[\([i, 
                j]\)] = \((∏\+\(m = 1\)\%NS\((KroneckerDelta[
                      Basis[\([j, m]\)] - 
                        Basis[\([i, m]\)]])\))\)*
			\((Sum[\((spin - 
                        Basis[\([j, l]\)])\)*\((spin - 
                        Basis[\([j, l + 1]\)])\), {l, 1, 
                    NS}])\)
]];\)\)

Off-diagonal = XY model

\(\(\(For[i = 1, 
        i < DimH + 1, \(i++\), 
For[j = 1, 
          j < DimH + 1, \(j++\), 
For[k = 1, 
            k < NS + 1, \(k++\), 
HilfsElement = 
              Table[Basis[\([j, n]\)], {n, 1, 
                  NS + 1}]; 
HilfsElement[\([k]\)] = 
              HilfsElement[\([k]\)] - 
                1; 
HilfsElement[\([k + 1]\)] = 
              HilfsElement[\([k + 1]\)] + 1; 
If[
              k == 1, 
              HilfsElement[\([NS + 1]\)] = 
                HilfsElement[\([1]\)]]; 
If[k == NS, 
              HilfsElement[\([1]\)] = 
                HilfsElement[\([NS + 1]\)]]; 
HMatrix[\([i, 
                  j]\)] = 
              HMatrix[\([i, j]\)] + 
                0.5*
Product[
                    KroneckerDelta[
                      Basis[\([i, m]\)] - HilfsElement[\([m]\)]], {m, 1, 
                      NS}]*
\@\(\((Basis[\([j, k]\)])\)*\((2*spin - Basis[\([j, k]\)] + 1)\)\)*
\@\(\((2*spin - Basis[\([j, k + 1]\)])\)*\((Basis[\([j, k + 1]\)] + 1)\)\)
]]];\)\(
\)
    \)\)
For[i = 1, 
        i < DimH + 1, i++, 
For[j = 1, 
          j < DimH + 1, j++, 
For[k = 1, 
            k < NS + 1, k++, 
HilfsElement = 
              Table[Basis[[j, n]], {n, 1, 
                  NS + 1}]; 
HilfsElement[[k]] = 
              HilfsElement[[k]] + 
                1; 
HilfsElement[[k + 1]] = 
              HilfsElement[[k + 1]] - 1; 
If[
              k == 1, 
              HilfsElement[[NS + 1]] = 
                HilfsElement[[1]]]; 
If[k == NS, 
              HilfsElement[[1]] = 
                HilfsElement[[NS + 1]]]; 
HMatrix[[i, 
                  j]] = 
              HMatrix[[i, j]] + 
                0.5*
Product[
                    KroneckerDelta[
                      Basis[[i, m]] - HilfsElement[[m]]], {m, 1, 
                      NS}]*
(Sqrt[(Basis[[j, k + 1]])*(2*spin - Basis[[j, k + 1]] + 1)])*
(Sqrt[(2*spin - Basis[[j, k]])*(Basis[[j, k]] + 1)])
]]];

If you are desperate, have a look at the full matrix!

Eigenvalues of H

EV = Sort[Chop[Eigenvalues[HMatrix]]]

Now, let's do the same in subspaces H(M)

M = 0;

MBasis = 
        Table[{0, 0, 0, 0, 0, 0, 0}, {n, 1, 
            1}];

For[i = 1, 
        i < DimH + 
            1, i++, 
ThisA = (∑_(j = 1)^NS) 
              Basis[[i, j]]; 
If[
          ThisA == NS*spin - M, 
          AppendTo[MBasis, 
            Basis[[i]]]]];

MBasis = Delete[MBasis, 1];

Check new basis

MatrixForm[MBasis]

DimHM = Length[MBasis]

Create Hamilton matrix in the subspace

Diagonal = Ising model

HMMatrix = Table[0.0, {i, 1, DimHM}, {j, 1, DimHM}];
\(\(For[i = 1, 
        i < DimHM + 1, \(i++\), 
For[j = 1, 
          j < DimHM + 
              1, \(j++\), 
		HMMatrix[\([i, 
                j]\)] = \((∏\+\(m = 1\)\%NS\((KroneckerDelta[
                      MBasis[\([j, m]\)] - 
                        MBasis[\([i, m]\)]])\))\)*
			\((Sum[\((spin - 
                        MBasis[\([j, l]\)])\)*\((spin - 
                        MBasis[\([j, l + 1]\)])\), {l, 1, 
                    NS}])\)
]];\)\)

Off-diagonal = XY model

\(\(For[i = 1, 
        i < DimHM + 1, \(i++\), 
For[j = 1, 
          j < DimHM + 1, \(j++\), 
For[k = 1, 
            k < NS + 1, \(k++\), 
HilfsElement = 
              Table[MBasis[\([j, n]\)], {n, 1, 
                  NS + 1}]; 
HilfsElement[\([k]\)] = 
              HilfsElement[\([k]\)] - 
                1; 
HilfsElement[\([k + 1]\)] = 
              HilfsElement[\([k + 1]\)] + 1; 
If[
              k == 1, 
              HilfsElement[\([NS + 1]\)] = 
                HilfsElement[\([1]\)]]; 
If[k == NS, 
              HilfsElement[\([1]\)] = 
                HilfsElement[\([NS + 
                      1]\)]]; 
HMMatrix[\([i, j]\)] = 
              HMMatrix[\([i, j]\)] + 
                0.5*
Product[
                    KroneckerDelta[
                      MBasis[\([i, m]\)] - HilfsElement[\([m]\)]], {m, 1, 
                      NS}]*
\@\(\((MBasis[\([j, k]\)])\)*\((2*spin - MBasis[\([j, k]\)] + 1)\)\)*
\@\(\((2*spin - MBasis[\([j, k + 1]\)])\)*\((MBasis[\([j, k + 1]\)] + 1)\)\)
]]];\)
\)
For[i = 1, 
        i < DimHM + 1, i++, 
For[j = 1, 
          j < DimHM + 1, j++, 
For[k = 1, 
            k < NS + 1, k++, 
HilfsElement = 
              Table[MBasis[[j, n]], {n, 1, 
                  NS + 1}]; 
HilfsElement[[k]] = 
              HilfsElement[[k]] + 
                1; 
HilfsElement[[k + 1]] = 
              HilfsElement[[k + 1]] - 1; 
If[
              k == 1, 
              HilfsElement[[NS + 1]] = 
                HilfsElement[[1]]]; 
If[k == NS, 
              HilfsElement[[1]] = 
                HilfsElement[[NS + 
                      1]]]; 
HMMatrix[[i, j]] = 
              HMMatrix[[i, j]] + 
                0.5*
Product[
                    KroneckerDelta[
                      MBasis[[i, m]] - HilfsElement[[m]]], {m, 1, 
                      NS}]*
(Sqrt[(MBasis[[j, k + 1]])*(2*spin - MBasis[[j, k + 1]] + 1)])*
(Sqrt[(2*spin - MBasis[[j, k]])*(MBasis[[j, k]] + 1)])
]]];

Eigenvalues of H in H(M)

EV = Sort[Chop[Eigenvalues[HMMatrix]]]