ED - Yan's group

For exact diagonalization (ED) studies, the foundational approach involves delving into several fundamental models for beginners. Examples include the lattice fermion Hubbard Model and quantum spin models such as the Ising Model and the XXZ Model (Heisenberg Model).

Our emphasis here is to introduce how to construct the Hamiltonian, commencing with the Heisenberg standard model and extending towards the XXZ model. The Ising Model assumes a more simplified form due to its absence of off-diagonal interaction terms.

1D Heisenberg chain

the Hamiltonian gives:
$$\begin{array}{}& \begin{array}{c}H=J\sum _{<i,j>}{\mathbf{S}}_i\cdot {\mathbf{S}}_j\\ \to H=J\sum _{⟨i,j⟩}(S_i^x{S}_j^x+S_i^y{S}_j^y+S_i^z{S}_j^z)\end{array}\end{array}$$

According to the relationships:$S^x=\frac{1}{2}(S^{+}+S^{-}),S^y=\frac{1}{2i}(S^{+}-S^{-})$

By expressing the Hamiltonian in terms of the $S^{+},S^{-}$:

$$H=[\frac{1}{2}(S_i^{+}{S}_j^{-}+S_i^{-}{S}_j^{+})+S_i^z{S}_j^z]$$

Given the presence of two sites representing spins, the size of the Hilbert space is determined by

$$D=2(\text{Two spin orientations.}{)}^{2(\text{Two sites})}=4$$,

Consequently, the Hamiltonian matrix for the system is a 4x4 matrix. Describing the two sites is straightforward, and we typically represent the Hamiltonian solely in operator form. From a matrix perspective, the interactions between operators are connected through tensor products to link matrices.

$$\text{matrix representing: }H=[\frac{1}{2}(S_1^{+}\otimes S_2^{-}+S_1^{-}\otimes S_2^{+})+S_1^z\otimes S_2^z]$$

The fundamental idea behind Exact Diagonalization (ED) is as follows: After obtaining the matrix representation of the Hamiltonian, numerical diagonalization is performed on the Hamiltonian matrix. In theory, this process allows us to obtain all the eigenvalues and eigenvectors associated with the Hamiltonian, providing a comprehensive description of the system.

The method we are about to introduce below is the simplest and most foundational, relying on intuitive understanding, without involving the treatment of larger system sizes.

Hamiltonian

Since we understand that the Hamiltonian matrix is in the form of the tensor product of operators, for a small number of sites, such as a 2-site matrix being 4x4, a 3-site matrix being 8x8, and a 4-site matrix being 16x16, the matrix elements can be manually calculated. Therefore, we can cross-verify the manually calculated matrix form with the matrix form provided by the code.

The manually calculated matrix form of the Hamiltonian is very easy to obtain and we will now focus on explaining the code implementation.

Firstly, we pre-store the matrix forms of the operators that constitute the Hamiltonian:
$$\begin{array}{}& \begin{array}{c}{S}^z=\frac{1}{2}{\sigma }^z=\frac{1}{2}\left(\begin{array}{cc}1& \\ & -1\end{array}\right),S^{+}=\left(\begin{array}{cc}0& 1\\ 0& 0\end{array}\right),I=\left(\begin{array}{cc}1& \\ & 1\end{array}\right),\hslash =1\end{array}\end{array}$$

sz = 0.5*[1.0,0;0,-1.0];
sp = [0,1;0,0];
id = eye(2,2);

sz =
0.5000         0
     0   -0.5000

sp =
 0     1
 0     0

id =
 1     0
 0     1

We do not redundantly store the operator $S^{-}$, as it is evident that its matrix form is merely the transpose of $S^{+}$. In MATLAB, adding a ‘ to a matrix, such as A’, denotes the transpose of the matrix A.

Before we construct the Hamiltonian matrix, it is crucial to create(the matrix) and initialize:

H = zeros(2,2)
Sp0 = sp
Sz0 = sz
Id = eye(2,2)

Sp1 = Sp0
Sz1 = Sz0

H =
     0     0
     0     0


Sp0 =
     0     1
     0     0


Sz0 =
    0.5000         0
         0   -0.5000


Id =
     1     0
     0     1


Sp1 =
     0     1
     0     0


Sz1 =
    0.5000         0
         0   -0.5000

Notice that we have additionally prepared operator matrices Sp0, Sp1, Sz0, Sz1, Id (MATLAB is case-sensitive, hence Id and id are distinct variables). These matrices have been initialized, and the form of the initialization is provided by the results.

With all the preparatory work completed, we can now proceed to the main loop for constructing the Hamiltonian.

Evidently, the core of this loop is:

for site=1:L-1
	...
end

The variable site, representing the site point coordinates, concludes the loop at L-1 instead of L. This is because, in practice, when dealing with site = i, then $j=site+1$. If the loop were to end at L, the value of j would exceed the length.

Now, we can contemplate how to create the code for our first 4x4 Hamiltonian matrix. It is:

$$H=[\frac{1}{2}(S_1^{+}\otimes S_2^{-}+S_1^{-}\otimes S_2^{+})+S_1^z\otimes S_2^z]$$

We can denote $S_1^{+}\otimes S_2^{-}$ as h. In MATLAB, the tensor product of two matrices can be represented using the kron(A,B) function. Therefore, $h=S_1^{+}\otimes S_2^{-}$. Similarly, it is evident that $h^T=(S_1^{+}\otimes S_2^{-}{)}^T={(S_1^{+})}^T\otimes {(S_2^{-})}^T=S_1^{-}\otimes S_2^{+}$.

h = 0.5*kron(Sp1, sp');
H = kron(H, id) + h + h' + kron(Sz1,sz);

Here, both Sp1 = sp and Sz1 = sz are 2x2 matrices. After the tensor product, they become 4x4 matrices, i.e., kron(H, id), h, h', and kron(Sz1, sz) are all 4x4 matrices. Matrix addition relies on having the same dimensions. The first term in the Hamiltonian matrix serves the purpose of expanding the initialized Hamiltonian matrix H from 2x2 to 4x4, maintaining consistency with the dimensions of the subsequent terms.

2-spin system

If we consider only 2-spin system, it is unsurprising that the Hamiltonian is immediately fully constructed, here we print the result:

H =
    0.2500         0         0         0
         0   -0.2500    0.5000         0
         0    0.5000   -0.2500         0
         0         0         0    0.2500

4*H =
     1     0     0     0
     0    -1     2     0
     0     2    -1     0
     0     0     0     1

Consistent with our manual derivation. The matrices below are four times the previous ones, which can be interpreted as factoring out J/4.

3-spin system

For 3-spin system, to maintain consistency in the dimensions of the Hamiltonian matrix, it should be:
$$\begin{array}{}& \begin{array}{c}H=[\frac{1}{2}(S_1^{+}\otimes S_2^{-}\otimes I_3+S_1^{-}\otimes S_2^{+}\otimes I_3)+S_1^z\otimes S_2^z\otimes I_3\\ +\frac{1}{2}(I_1\otimes S_2^{+}\otimes S_3^{-}+I_1\otimes S_2^{-}\otimes S_3^{+})+I_1\otimes S_2^z\otimes S_3^z\\ +\frac{1}{2}(S_1^{+}\otimes I_2\otimes S_3^{-}+S_1^{-}\otimes I_2\otimes S_3^{+})+S_1^z\otimes I_2\otimes S_3^z]\end{array}\end{array}$$

We call it Eq.7. It is not difficult to comprehend. The first row signifies the action of the interaction operators $S_i^{\pm ,z}{S}_j^{\mp ,z}$ on lattice points where i=1 and j=2. However, to maintain dimensional consistency, it is as if the identity operator acts on the third lattice point and contributes to the matrix. (One can contemplate the consequences if we express H solely as $[\frac{1}{2}(S_1^{+}\otimes S_2^{-}+S_1^{-}\otimes S_2^{+})+S_1^z\otimes S_2^z]+[\frac{1}{2}(S_2^{+}\otimes S_3^{-}+S_2^{-}\otimes S_3^{+})+S_2^z\otimes S_3^z]+..$
There would arise significant issues.) The Hamiltonian matrix is fundamentally expressed based on the chosen basis, and therefore, the order of the basis for each row and column has implications.)

It is noteworthy that the last row indicates the action of two-body interactions on points i=3 and j=1, suggesting the adoption of periodic boundary conditions in this scenario.

Now, let’s delve into each row specifically:

$$\begin{array}{}& \begin{array}{c}{H}_1=\frac{1}{2}(S_1^{+}\otimes S_2^{-}\otimes I_3+S_1^{-}\otimes S_2^{+}\otimes I_3)+S_1^z\otimes S_2^z\otimes I_3\\ \to H_1=H_{old}\otimes I\end{array}\end{array}$$

we realize that the Hamiltonian $H_1$ in the first row, factoring out $\otimes I_3$, is essentially equivalent to the Hamiltonian for two sites that we previously calculated, denoted as $H_{old}=H_{2-spin}$. Consequently, in the code:

% when site=1 for 2-spin
	h = 0.5*kron(Sp1, sp');
	H = kron(H, id) + h + h' + kron(Sz1,sz);

% for 3-spin
	H_1 = kron(H,id)

We have implemented the first row of the total Hamiltonian(Eq.7).

Now, the second row of Hamiltonian(Eq.7):

$$H_2=\frac{1}{2}(I_1\otimes S_2^{+}\otimes S_3^{-}+I_1\otimes S_2^{-}\otimes S_3^{+})+I_1\otimes S_2^z\otimes S_3^z$$

Sp1 = kron(Id, sp); 
Sz1 = kron(Id, sz); 

% outside:
h = 0.5*kron(Sp1, sp');
H_2 = h + h' + kron(Sz1,sz);

Still remember that our Sp1 and Sz1 were initialized as sp and sz 2x2 matrices? After the operation Sp1 = kron(Id, sp), it represents: $I_1\otimes S_2^{+}$. Similarly, Sz1 = kron(Id, sz) represents $I_1\otimes S_2^z$. Considering them as a whole, the term h = 0.5*kron(Sp1, sp') now represents $1/2(I_1\otimes S_2^{+}\otimes S_3^{-})$, and h' represents $1/2(I_1\otimes S_2^{-}\otimes S_3^{+})$. Finally, the term kron(Sz1,sz) represents $I_1\otimes S_2^z\otimes S_3^z$. This way, we have completed the construction of the second row of the Hamiltonian.

For the third row of Hamiltonian(Eq.7) , we handle the last point, where i = 3:

Sp0 = kron(Sp0, id);
Sz0 = kron(Sz0, id);

% outside:
h = 0.5*kron(Sp0, sp'); 
H = H + h + h' + kron(Sz0,sz); 

$$H_3=\frac{1}{2}(S_1^{+}\otimes I_2\otimes S_3^{-}+S_1^{-}\otimes I_2\otimes S_3^{+})+S_1^z\otimes I_2\otimes S_3^z$$

Here, Sp0 = kron(Sp0, id) with the initialization Sp0 = sp represents $S_1^{+}\otimes I_2$. Therefore, h = 0.5*kron(Sp0, sp') corresponds to $1/2(S_1^{+}\otimes I_2\otimes S_3^{-})$. Similarly, $S_1^z\otimes I_2\otimes S_3^z$ represents the treatment for the last point i=3. This completes the construction of the three components of the Hamiltonian for the 3-spin system under periodic boundary conditions.

Below is the complete code, composed of the snippets above, with some logical changes for functionality linkage. With the presence of comments, it should be comprehensible for readers:

% 3-spin example:
for site=1:L-1 
	
    h = 0.5*kron(Sp1, sp') 
    H = kron(H, id) + h + h' + kron(Sz1,sz)
    % If L = 2, the process concludes at this point. The code above will only execute once.
    % If L = 3, The first term H = kron(H, id) represents H_1 (Regardless of the number of points (L > 2), it always represents the Hamiltonian for the first row H_1 = H_old otimes I.)
    % ..h + h' + kron(Sz1,sz) represents H_2
    
    
    
		% Here is the code for extending the matrix when L > 2.
    if L>2
        % The following is our handling process for the last point: i = L-1, j = L = 1
        if site == L-1
            h = 0.5*kron(Sp0, sp');
            H = H + h + h' + kron(Sz0,sz); % H_3
        end

			% When L > 2 but site = 1, we handle matrices before entering the next iteration:site=L-1=2
        if site < L-1
            Sp1 = kron(Id, sp); 
            Sz1 = kron(Id, sz);

            Sp0 = kron(Sp0, id);
            Sz0 = kron(Sz0, id);
            Id = eye(size(Sp0)); 
        end
    end
end

Finally, let’s output the Hamiltonian matrix and compare it with the manually derived form.

H =
    0.7500         0         0         0         0         0         0         0
         0   -0.2500    0.5000         0    0.5000         0         0         0
         0    0.5000   -0.2500         0    0.5000         0         0         0
         0         0         0   -0.2500         0    0.5000    0.5000         0
         0    0.5000    0.5000         0   -0.2500         0         0         0
         0         0         0    0.5000         0   -0.2500    0.5000         0
         0         0         0    0.5000         0    0.5000   -0.2500         0
         0         0         0         0         0         0         0    0.7500
         
4*H =
     3     0     0     0     0     0     0     0
     0    -1     2     0     2     0     0     0
     0     2    -1     0     2     0     0     0
     0     0     0    -1     0     2     2     0
     0     2     2     0    -1     0     0     0
     0     0     0     2     0    -1     2     0
     0     0     0     2     0     2    -1     0
     0     0     0     0     0     0     0     3

---

The code demonstrated above is also applicable for L = 4, L = 5, L = 6, … since the implementation logic remains the same and covers cases with different numbers of lattice sites. I won’t go into excessive details here. Instead, let’s just print out the Hamiltonian matrix for a 4-spin system. Readers can compare it with the manually derived form:

$16\times 16$ Hamiltonian matrix of standard Heisenberg Model:

H = J/4 x 

     4     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     2     0     0     0     0     0     2     0     0     0     0     0     0     0
     0     2     0     0     2     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     2     0     0     0     0     2     0     0     0     0     0
     0     0     2     0     0     0     0     0     2     0     0     0     0     0     0     0
     0     0     0     2     0    -4     2     0     0     2     0     0     2     0     0     0
     0     0     0     0     0     2     0     0     0     0     2     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0     0     2     0     0     2     0
     0     2     0     0     2     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     2     0     0     0     0     2     0     0     0     0     0
     0     0     0     2     0     0     2     0     0     2    -4     0     2     0     0     0
     0     0     0     0     0     0     0     2     0     0     0     0     0     2     0     0
     0     0     0     0     0     2     0     0     0     0     2     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0     0     2     0     0     2     0
     0     0     0     0     0     0     0     2     0     0     0     0     0     2     0     0
     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     4

Appendix A: tensor products rules

For more details, please refer to Wikipedia.

Tensor product does not require matrices to have the same dimensions. Suppose there is a 2x2 matrix A and a 2x2 matrix B for this demonstration.

$$\begin{array}{}& \begin{array}{c}A=\left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right],B=\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]\end{array}\end{array}$$

$$\begin{array}{}& \begin{array}{c}A\otimes B=\left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right]\otimes \left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]\end{array}\end{array}$$
$$\begin{array}{}& \begin{array}{c}=\left[\begin{array}{cc}1\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]& 0\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]\\ 0\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]& -1\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]\end{array}\right]=\left[\begin{array}{cccc}& 1& & \\ 1& & & \\ & & & -1\\ & & -1& \end{array}\right]\end{array}\end{array}$$
$$\begin{array}{}& \begin{array}{c}B\otimes A=\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]\otimes \left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right]\end{array}\end{array}$$
$$\begin{array}{}& \begin{array}{c}=\left[\begin{array}{cc}0\left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right]& 1\left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right]\\ 1\left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right]& 0\left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right]\end{array}\right]=\left[\begin{array}{cccc}& & 1& \\ & & & -1\\ 1& & & \\ & -1& & \end{array}\right]\end{array}\end{array}$$

Appendix B: Pauli matrices

Pauli matrix:
$$\begin{array}{}& \begin{array}{c}{\sigma }^z=\left(\begin{array}{cc}1& \\ & -1\end{array}\right),{\sigma }^x=\left(\begin{array}{cc}& 1\\ 1& \end{array}\right),{\sigma }^y=\left(\begin{array}{cc}& -i\\ i& \end{array}\right)\end{array}\end{array}$$

The relationship between $S=1/2$ operators and Pauli matrices is as follows:

$$\overrightarrow{S}=\frac{\hslash }{2}\overrightarrow{\sigma },\overrightarrow{\sigma }=({\sigma }^x,{\sigma }^y,{\sigma }^z)$$

We typically assume $\hslash =1$, and given the relationship between the spin operator with $S=1/2$ and the Pauli matrices, it is possible to derive the matrix forms of the raising and lowering spin operators based on these relationships.

$$\begin{array}{}& \begin{array}{c}{S}^{+}=\frac{\hslash }{2}(\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)+i\left(\begin{array}{cc}0& -i\\ i& 0\end{array}\right))=\frac{\hslash }{2}\left(\begin{array}{cc}0& 1-1\\ 1-i& 0\end{array}\right)=\hslash \left(\begin{array}{cc}0& 1\\ 0& 0\end{array}\right)\end{array}\end{array}$$

$$\begin{array}{}& \begin{array}{c}{S}^{-}=\frac{\hslash }{2}(\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)-i\left(\begin{array}{cc}0& -i\\ i& 0\end{array}\right))=\frac{\hslash }{2}\left(\begin{array}{cc}0& 1-1\\ 1+1& 0\end{array}\right)=\hslash \left(\begin{array}{cc}0& 0\\ 1& 0\end{array}\right)\end{array}\end{array}$$

Clearly,

$$\begin{array}{}& \begin{array}{c}{S}^{+}|\uparrow ⟩=0,S^{+}|\downarrow ⟩=|\uparrow ⟩\\ S^{-}|\uparrow ⟩=|\downarrow ⟩,S^{-}|\downarrow ⟩=0\end{array}\end{array}$$

Appendix C: external field

1D Longitudinal field

sh = zeros(2,2);
idnew = id;
szh = sz;
for site=1:L-1 %L-1

    sh = kron(sh,id) + hb * kron(szh,id);

    if L>2
        if site == L-1
            sh = sh + hb*kron(idnew,sz);
        end


        if site < L-1
            szh = kron(idnew, sz); 
            idnew = eye(size(szh)); 
        end
    else %(For L=2)
     		sh = sh + hb * kron(id,sz);
    end

end

1D Transverse field

sx = [0,1.0;1.0,0];%notice here is Pauli Matrix
iidnew = id;
sxg = sx;
sgamma = zeros(2,2);
for site=1:L-1 %L-1

    sgamma = kron(sgamma,id) + gamma * kron(sxg,id);

    if L>2

        if site == L-1
            sgamma = sgamma + gamma * kron(iidnew,sx);
        end


        if site < L-1
            sxg = kron(iidnew, sx); 
            iidnew = eye(size(sxg)); 
        end
    else 
    sgamma = sgamma + gamma * kron(id,sx);
    end

end