Tutorial: 12-site Kitaev cluster¶
Here we will work through an example of using the library to calculate the ground state properties of the spin-1 Kitaev Hamiltonian with single-ion anisotropy on a lattice with \(L=12\) sites. The full Hamiltonian for this model is given by:
Here \(K_x\), \(K_y\), and \(K_z\) are the Kitaev couplings along the \(x\), \(y\), and \(z\) bond directions of the honeycomb lattice. The Hamiltonian includes a term describing single-ion anisotropy in the \([1, 1, 1]\) direction, with magnitude \(D\), which preserves the symmetry between \(x\), \(y\), and \(z\) bonds. For the spin-1 case, the Hilbert space has dimension \(3^L\), where \(L\) is the total number of sites. However, using lattice symmetries we can greatly reduce the dimensions of the Hamiltonian (by a factor of \(L\)), and through efficient construction and storage of the sparse matrices \(H_{x}\), \(H_{y}\), \(H_{z}\), and \(H_{D}\), it becomes feasible to study clusters with \(L \geq 18\) sites.
1. Importing modules¶
Let’s start by importing the modules we will need. After specifying the size of our lattice, we will be using the symmetry functions
get_neighbors, lattice_translations and get_representative_states. These find all nearest-neighbor pairs, the configurations
of the lattice which are equivalent via symmetry, and the set of representative states we need to keep (which greatly reduces the size
of the Hilbert space and therefore the dimensions of our Hamiltonian matrix).
Then we will use a set of functions to efficiently construct and store the components of the Hamiltonian: construct_Hxx, construct_Hyy, construct_Hzz, construct_HD,
and finally construct_hamiltonian to reproduce the full Hamiltonian. Then we will obtain the ground state and calculate its energy per site \(E_0\),
bipartite entanglement entropy, and the expectation value of the squared spin operator \((S^x + S^y + S^z)^2\). For this we will need to import the functions
ground_state, s_squared, and entanglement_entropy:
import numpy as np
import scipy
import gc
from kitaev_clusters.symmetry_functions import get_neighbors, lattice_translations, get_representative_states
from kitaev_clusters.hamiltonian_functions import construct_Hxx, construct_Hyy, construct_Hzz, construct_HD, construct_hamiltonian
from kitaev_clusters.ground_state_functions import ground_state, s_squared, entanglement_entropy
2. Defining lattice parameters¶
Each unit cell of the honeycomb lattice consists of two sites. For a lattice with \(L=12\) sites, we can choose a geometry with \(L_x=3\) unit cells
in the x-direction, and \(L_y=2\) unit cells in the y-direction. Then we will call the function get_neighbors to get a list of site pairs along each
of the three bond directions of the honeycomb lattice. After this, we call the function lattice_translations to get a list of all permutations of \(L\)
sites which are equivalent via honeycomb lattice symmetries:
L = 12 # Total number of lattice sites (L = Lx * Ly * 2)
Lx = 3 # Unit cells in horizontal direction
Ly = 2 # Unit cells in vertical direction
x_neighbors = get_neighbors(Lx, Ly, bond='x')
y_neighbors = get_neighbors(Lx, Ly, bond='y')
z_neighbors = get_neighbors(Lx, Ly, bond='z')
translation_order_list = lattice_translations(Lx, Ly)
3. Finding representative states¶
Here we will find the representative states we need to keep using the function get_representative_states, and construct and store the following arrays.
We find kept_ints, an array of the base-10 integers corresponding to the representative states. We also get state_map, an array containing the representative
state (in base-10) for each of the \(3^L\) possible states, i.e. a mapping from a state to its representative state. Each representative state has a number of states
equivalent to it by symmetry, however only a certain number of these will be unique. We store n_unique_list, which is an array containing the number of unique
mirror states for each representative state. Finally, the function also returns ndim, which is the total number of representative states we keep, i.e. the dimension
of the reduced Hilbert space.
kept_ints, state_map, n_unique_list, ndim = get_representative_states(L, translation_order_list,
ints_filename='kept_ints_L12',
states_filename='state_map_L12',
nunique_filename='n_unique_list_L12')
4. Loading representative states¶
This step is only necessary if you want to load from file the arrays of representative states, state map, and number of unique mirror states. However, for a given
lattice size, these arrays will be fixed and therefore only need to be calculated once. It is therefore recommended to load these arrays from file. Note that the
function get_representative_states will save these arrays in .npz format:
kept_ints = np.load('kept_ints_L12.npy')
state_map = np.load('state_map_L12.npy')
n_unique_list = np.load('n_unique_list_L12.npy')
5. Constructing components of Hamiltonian¶
Here we pass our previously loaded arrays to the functions construct_Hxx, construct_Hyy, construct_Hzz, and construct_HD to efficiently construct
the components of the Hamiltonian as sparse matrices in List of List (LIL) format. These matrices are then converted to Compressed Sparse Row (CSR) format
for reduced memory requirements (although this can be disabled by setting as_csr=False), and saved in .npz format. We can specify the filename for each
matrix to be saved.
Hxx = construct_Hxx(L, x_neighbors, kept_ints, state_map, n_unique_list, filename='Hxx_L12', as_csr=True)
Hyy = construct_Hyy(L, y_neighbors, kept_ints, state_map, n_unique_list, filename='Hyy_L12', as_csr=True)
Hzz = construct_Hzz(L, z_neighbors, kept_ints, filename='Hzz_L12', as_csr=True)
HD = construct_HD(L, kept_ints, state_map, n_unique_list, filename='HD_L12', as_csr=True)
6. Loading components of Hamiltonian¶
For a given lattice size, the matrices \(H_x\), \(H_y\), \(H_z\), and \(H_D\) are fixed. It is therefore recommended to only construct them once then load from file when needed. Since these are sparse matrices saved in .npz format, we can load them as follows:
Hxx = scipy.sparse.load_npz('Hxx_L12.npz')
Hyy = scipy.sparse.load_npz('Hyy_L12.npz')
Hzz = scipy.sparse.load_npz('Hzz_L12.npz')
HD = scipy.sparse.load_npz('HD_L12.npz')
However, it is not necessary to load each component of the Hamiltonian separately if you just want to use them to construct the full Hamiltonian. That
is performed using the function construct_hamiltonian as shown in the next step.
7. Specifying couplings, anisotropy, and constructing full Hamiltonian¶
Now we can specify values for the Kitaev couplings \(K_x\), \(K_y\), and \(K_z\), as well as the anisotropy parameter \(D\):
Kx = 1 # Kitaev coupling along x-bonds
Ky = 1 # Kitaev coupling along y-bonds
Kz = 1 # Kitaev coupling along z-bonds
D = 0.0 # [1, 1, 1] single-ion anisotropy
The function construct_hamiltonian can now be called, passing these parameters as arguments, to return the full Hamiltonian matrix \(H\).
This function loads each of the components of the Hamiltonian from file, with the filenames of these .npz files also passed as arguments:
H = construct_hamiltonian(Kx, Ky, Kz, D, 'Hxx_L12.npz', 'Hyy_L12.npz', 'Hzz_L12.npz', 'HD_L12.npz')
However, if you are not loading the components of the Hamiltonian from file, for example if you already have constructed the matrices \(H_x\), \(H_y\), \(H_z\), and \(H_D\) in your workspace or notebook, simply create the full Hamiltonian matrix as follows:
H = (Kx * Hxx) + (Ky * Hyy) + (Kz * Hzz) + (D * HD)
8. Finding ground sate and physical properties¶
Now we can pass the full Hamiltonian matrix \(H\) as an argument to the function ground_state to obtain the ground state psi_0 and
the ground state energy per site E_0. We can also pass the ground state psi_0 as an argument to the functions s_squared and
entanglement_entropy to get the expectation value of \((S^x + S^y + S^z)^2\), and the bipartite entanglement entropy:
E_0, psi_0 = ground_state(H, L)
Ss = s_squared(HD, L, psi_0)
entropy = entanglement_entropy(L, psi_0, kept_ints, state_map, n_unique_list)
Saving these results to a .txt file is also straightforward:
f = open('results.txt', 'w')
f.write('D_111, E_gs per site, <(Sx + Sy + Sz)^2>, entropy \n')
f.write(f'{D}, {E_0}, {Ss}, {entropy} \n')
f.close()
If desired, you may want to saved the ground state psi_0 after calling the function ground_state, and later load it when calculating
the entanglement entropy or \(\langle(S^x + S^y + S^z)^2\rangle\). This can be performed with the commands:
# Saving ground state as .npy file:
np.save('psi_0', psi_0)
# Loading ground state:
np.load('psi_0.npy')
9. Freeing memory after calculations complete¶
If desired, you may want to free up memory after your calculations are complete. If you no longer need to store the matrices you have loaded or the ground state, this can be done as follows:
del Hxx
del Hyy
del Hzz
del HD
del H
del psi_0
gc.collect()