Eigen-decomposition

Let's construct the hopping term \(\sum_{\sigma } {\hat{f} }_{2\sigma }^{\dagger } {\hat{f} }_{1\sigma } +{\hat{f} }_{1\sigma }^{\dagger } {\hat{f} }_{2\sigma }\) acting on two spinful fermionic sites.

% for site s00
F1 = F; F1.info.itags = { 's00', 's00*', 'op*'};
E1 = I.E; E1.info.itags = { 's00', 's00*'};
% for site s01
F2 = F; F2.info.itags = { 's01', 's01*', 'op*'};
E2 = I.E; E2.info.itags = { 's01', 's01*'};
Z2 = Z; Z2.info.itags = { 's01', 's01*'};
A = getIdentity(E1,2,E2,2, 'A01*',[1 3 2]);

H = contract(A, '!2*',{F1, '!1',{F2, '!2*',{Z2, '!1',A}}}) + ...
    contract(A, '!2*',{F1, '!2*',{Z2, '!1',{F2, '!1',A}}}) + ...
    getIdentity(A,2) * 1e-40

The first line of defining H means \(\sum_{\sigma } {\hat{f} }_{2\sigma }^{\dagger } {\hat{f} }_{1\sigma }\), and the second line means its Hermitian conjugate. And in the third line, we added the identity multiplied by a small number, to let H have all the sectors (that amount to 16 dimensional space).

celldisp(H.data)

ans{1} =
   1.0000e-40


ans{2} =
   0.0000    1.0000
   1.0000    0.0000

ans{3} =
    0.0000   -1.4142          0
   -1.4142    0.0000    -1.4142
         0   -1.4142     0.0000

ans{4} =
   1.0000e-40


ans{5} =
    0.0000   -1.0000
   -1.0000    0.0000

ans{6} =
   1.0000e-40

The eigenvalues and eigenvectors of H can be obtained by eig which is the wrap-up of eigQS.

[V,D] = eig(H)

V =
    Q: 6x [2 2] having 'A,SU2',      { A01, A01* }   
 data: 2-D double (784 bytes)       10 x 10 => 16 x 16

    1. 1x1        | 1x1       [ -2 0 ; -2 0 ]          1.
    2. 2x2        | 2x2       [ -1 1 ; -1 1 ]  32 B      {1.414}
    3. 3x3        | 1x1       [  0 0 ;  0 0 ]  72 B          
    4. 1x1        | 3x3       [  0 2 ;  0 2 ]          1.  {1.732}
    5. 2x2        | 2x2       [  1 1 ;  1 1 ]  32 B        {1.414}
    6. 1x1        | 1x1       [  2 0 ;  2 0 ]          1.

D =
    Q: 6x [2 2] having 'A,SU2',      { A01, A01* }   
 data: 2-D double (704 bytes)       6 x 10 => 10 x 16

    1. 1x1        | 1x2       [ -2 0 ; -2 0 ]          1e-40
    2. 2x2        | 2x2       [ -1 1 ; -1 1 ]  16 B      {1.414}
    3. 1x3        | 1x1       [  0 0 ;  0 0 ]  24 B          
    4. 1x1        | 3x3       [  0 2 ;  0 2 ]          1e-40  {1.732}
    5. 1x2        | 2x2       [  1 1 ;  1 1 ]  16 B        {1.414}
    6. 1x1        | 1x1       [  2 0 ;  2 0 ]          1e-40

F and D are QSpace objects. Each data of V is the unitary matrix whose columns are eigenvectors:

celldisp(V.data)

ans{1} =
   1


ans{2} =
  -0.7071    0.7071
   0.7071    0.7071

ans{3} =
   -0.5000   -0.7071    -0.5000
   -0.7071   -0.0000     0.7071
   -0.5000    0.7071    -0.5000

ans{4} =
   1


ans{5} =
   -0.7071   -0.7071
   -0.7071   -0.7071

ans{6} =
   1

We can check the unitarity of V by

V2 = contract(V, '!2*',V);
sameas(V2,getIdentity(A,2))

ans = logical
   1

Each data of D is the row vector of eigenvalues:

celldisp(D.data)

ans{1} =
   1.0000e-40


ans{2} =
  -1.0000    -1.0000


ans{3} =
  -2.0000    -0.0000     2.0000


ans{4} =
   1.0000e-40


ans{5} =
  -1.0000    -1.0000


ans{6} =
   1.0000e-40

(Quick exercise: Explain the eigenvalues.)

To make D as an operator representing a diagonal matrix,

D2 = diag(D)

V =
    Q: 6x [2 2] having 'A,SU2',      { A01, A01* }   
 data: 2-D double (784 bytes)       10 x 10 => 16 x 16

    1. 1x1        | 1x1       [ -2 0 ; -2 0 ]          1e-40
    2. 2x2        | 2x2       [ -1 1 ; -1 1 ]  32 B      {1.414}
    3. 3x3        | 1x1       [  0 0 ;  0 0 ]  72 B          
    4. 1x1        | 3x3       [  0 2 ;  0 2 ]          1e-40  {1.732}
    5. 2x2        | 2x2       [  1 1 ;  1 1 ]  32 B        {1.414}
    6. 1x1        | 1x1       [  2 0 ;  2 0 ]          1e-40

celldisp(D2.data)

ans{1} =
   1.0000e-40


ans{2} =
  -1.0000         0
        0    1.0000

ans{3} =
   -2.0000         0          0
         0   -0.0000          0
         0         0     2.0000

ans{4} =
   1.0000e-40


ans{5} =
   -1.0000         0
         0   -1.0000

ans{6} =
   1.0000e-40

One may also directly use the original MEX function eigQS.

[E,Ieig] = eigQS(H);

Note that the syntax is a bit different from the wrap-up eig. E is two-column matrix whose first column is the energy eigenvalues (sorted in ascending order) and second column indicates the multiplet dimensions (i.e., degeneracy due to non-Abelian symmetry) associated with the eigenvalues. When only the Abelian symmetries are used, E becomes a column vector, without having the second column for the multiplet dimensions.

E

ans = 10x2
     -2.0000     1.0000
     -1.0000     2.0000
     -1.0000     2.0000
     -0.0000     1.0000
      0.0000     1.0000
      0.0000     1.0000
      0.0000     3.0000
      1.0000     2.0000
      1.0000     2.0000
      2.0000     1.0000

And Inrg is the struct that contains more result of the eigendecomposition, including the eigenvectors (.AK and .AT) and eigenvalues (.EK and .ET).

Ieig

ans = struct with fields:
       AK: [1x1 struct]
       AD: [1x1 struct]
       EK: [1x1 struct]
       ED: [1x1 struct]
       DB: [6x2 struct]
       NK: 10
   Etrunc: 0

Becasue of the MATLAB policy, the direct result of MEX functions should be of MATLAB built-in types, while the QSpace is the user-defined data type. Here .AK, .AT, .EK, and .ET are struct variables that are compatible with QSpace. So we wrap them up as QSpace objects:

Ieig.EK = QSpace(Ieig.EK);
Ieig.AK = QSpace(Ieig.AK);
Ieig.EK

ans =
    Q: 6x [2 2] having 'A,SU2',      { A01, A01* }   
 data: 2-D double (784 bytes)       10 x 10 => 16 x 16

    1. 1x1        | 1x1       [ -2 0 ; -2 0 ]          1e-40
    2. 2x2        | 2x2       [ -1 1 ; -1 1 ]  32 B      {1.414}
    3. 3x3        | 1x1       [  0 0 ;  0 0 ]  72 B          
    4. 1x1        | 3x3       [  0 2 ;  0 2 ]          1e-40  {1.732}
    5. 2x2        | 2x2       [  1 1 ;  1 1 ]  32 B        {1.414}
    6. 1x1        | 1x1       [  2 0 ;  2 0 ]          1e-40

Ieig.AK

ans =
    Q: 6x [2 2] having 'A,SU2',      { A01, A01* }   
 data: 2-D double (784 bytes)       10 x 10 => 16 x 16

    1. 1x1        | 1x1       [ -2 0 ; -2 0 ]          1.
    2. 2x2        | 2x2       [ -1 1 ; -1 1 ]  32 B      {1.414}
    3. 3x3        | 1x1       [  0 0 ;  0 0 ]  72 B          
    4. 1x1        | 3x3       [  0 2 ;  0 2 ]          1.  {1.732}
    5. 2x2        | 2x2       [  1 1 ;  1 1 ]  32 B        {1.414}
    6. 1x1        | 1x1       [  2 0 ;  2 0 ]          1.

We can set several options for eigQS, such as Nkeep (number of multiplets to be kept) and Etrunc (threshold energy such that the energy eigenvalues below the value are to be kept). For details, type: eigQS -?