I have a large matrix for a 1D zigzag edge model of an otherwise $3\times 3$ tight-binding Hamiltonian (3 basis functions, each corresponding to an atomic orbital), involving the variable $k_x$. The large matrix contains blocks of this $3\times 3$ Hamiltonian for each unit cell involved (I chose 100 unit cells). For different values of $k_x$, I numerically calculate eigenvalues and eigenvectors, and use the eigenvalues to plot the band structure shown below:
How do I identify which of the three orbitals the edge states are coming from, by using the numerical wavefunctions or otherwise?
What I've done so far: I used the following (HIGHLY INEFFICIENT) Python code to identify the indices of eigenvalues where the bands approximately cross. I was hoping to use this to find the eigenvectors corresponding to the point, and then identify which orbitals the edge states are coming from. My rationale is that the probability will suddenly jump to something else (or switch between bands) when the bands cross. I guess my biggest confusion is that I do not know which eigenvectors to choose, nor what they mean. That is, if my original matrix is tridiagonal, with $3\times 3$ blocks making up the system, what bases do the resulting eigenvectors have? In the case of energy eigenstates in just momentum space, I understand that each eigenvector corresponds to an energy wavefunction. But now I have the added complication of several unit cells.
find_evec_with_energy = 1; # crossing occurs around here
tol = 0.02
rows,cols = np.asarray(energies).shape # for each iteration of kx, I append the eigenvalues into the list called 'energies'
indices = []
for j in range(rows):
for k in range(cols):
if np.abs(np.asarray(energies)[j,k]-find_evec_with_energy)<tol:
indices.append(f"{j},{k}")
print("row,column of eigenvalue with desired energy")
print(indices)
Then, I guess I take choose the eigenvector corresponding to one of these indices, and then plot each triplet of rows, squared? I'm not sure.
n = 75 # index of unit cell where band crossing occurs??
wavefunction = evecs[:,n]
up,middle,down = wavefunction[::3],wavefunction[1::3],wavefunction[2::3]
plt.plot(np.abs(up)**2)
plt.show()
Sorry if this is a basic question, but I am only now starting to learn real-space material systems. Thanks. P.S. I'm not sure what tags to use.
EDIT: updated code:
First, loop over various momenta and get eigenvectors and eigenvalues:
energies = []
evecs = []
W = 100
for i in range(v_segs):
ham = H(vary_range[i],W)
ev,evc = np.linalg.eigh(ham)
energies.append(ev)
evecs.append(evc)
Then, find indices of energies near band crossing:
find_evec_with_energy = 1; # crossing occurs around here
tol = 0.02
rows,cols = np.array(energies).shape
print(rows)
print(cols)
indices = []
for j in range(rows):
for k in range(cols):
if np.abs(np.array(energies)[j,k]-find_evec_with_energy)<tol:
indices.append([j,k])
print("row,column of eigenvalue with desired energy")
print(indices) # outputs [[12, 100], [87, 100]]
Finally, try to find orbitals involved in band crossing:
def find_orbitals(index): #n is the row we want
ik,iband = index
wavefunction = np.array(evecs)[ik,:,iband]
w = wavefunction.conj() * wavefunction
up,middle,down = w[::3],w[1::3],w[2::3]
print(f"up,middle,down at index {index}")
tol = 0.01
if np.allclose(middle, 0, atol=tol):
print("edge states from up,down")
else:
print("not middle")
if np.allclose(up, 0, atol=tol):
print("edge states from middle,down")
else:
print("not up")
if np.allclose(down, 0, atol=tol):
print("edge states from up,middle")
else:
print("not down")
plt.plot(np.abs(up)**2)
plt.show()
plt.plot(np.abs(mid)**2)
plt.show()
plt.plot(np.abs(down)**2)
plt.show()
#list_of_indices = indices#[[12,100],[87,100]]
list_of_indices = [[12,100],[87,100]]
for q in range(len(list_of_indices)):
find_orbitals(list_of_indices[q])