# How to use wavefunctions/density to determine which orbitals lead to edge states?

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])


One way of determining this is using the projected density of states (P-DOS)

This resolves the DOS into specific orbitals thereby allowing you to discretize each orbitals weight for a specific energy.

$$\mathrm{PDOS}_\nu(E) = \sum_i \psi^*_{i,\nu} [\mathbf S | \psi_{i}\rangle]_\nu D(E-\epsilon_i)$$

Note here that $$|\psi_i\rangle$$ is the $$i$$th eigenvector and $$\nu$$ is the orbital index. The overlap matrix $$\mathbf S$$ is here for generality but for orthogonal basis it is the identity matrix, and finally $$D(E)$$ is the distribution function, Gaussian, Lorentzian etc.

The following analysis depends on what you want to figure out. If you are only interested in the weight of a given orbital for a given eigenstate, then doing $$\psi^*_{i,\nu}|\mathbf S|\psi_i\rangle$$ would give you the vector that gives you the weight for each orbital.
If you want this across the entire Brillouin zone you have to integrate over that, $$\sum_k w_k \mathrm{PDOS}_\nu(E)$$.

So here it seems your Python code is almost complete. You'll just want to sum the contribution for each block orbitals to get the full picture. Generally one may then plot the so-called fatbands which is basically the PDOS on top of the band-structure.

Note that the energy dependence on the P-DOS is not necessary if you are not interested in energy-resolved quantities.

If you want to track the eigenvectors in the Brillouin zone I would suggest you to sort according to the eigenvectors that resemble each other the most (besides a common phase-factor). I.e. something like $$\langle \psi_j|\mathbf S|\psi_i\rangle e^{i\theta}$$ where $$\theta$$ should be optimized to maximize the overlap between the states.

PS. It isn't totally clear to me what exactly you want to describe. Whether the PDOS is enough or whether you want some band-unfolding procedure isn't totally clear to me. Also, are you using the same small 3-orbital Hamiltonian repeated 100 times along one direction? This wouldn't be the best method.

EDIT: I have understood your question like this, you want to figure out how the edge states are dispersed across the orbitals in the system. And possibly which of the 3 orbitals (from the primary cell) have the largest contribution to the edge states.
I would suggest you start with a smaller system, say a width = 8 or less system such that the bands are easier to distinguish. Using a width of 100 only makes the spaghetti band-structure more difficult.
Since you already know which state is the edge-state I suggest you do as you first suggested. Have a specific index for the edge-state, then at some band-crossing (happening at a specific $$k$$) and do your PDOS like this. No need to complicate things if not necessary. :)

• Thank you for the thorough answer, +1. Would you mind clarifying "You'll just want to sum the contribution for each block orbitals to get the full picture."? For 100 unit cells, I get energies of size (100,300) and evecs of size (300,300). Let's say I get [12, 100] and [87, 100] as the [row,col] indices of the required energies. However, do I do wavefunction = evecs[:,100]? Or, for each orbital (of total 3) take two terms in the PDOS sum (= evecs[12,:] and = evecs[87,:])? I just want to find which two orbitals the edge states are coming from. Jun 22 at 15:46
• So far, the most promising has been to use wavefunction=evecs[:,n] for each row index n, to plot each orbital separately. Then, I see that one of these orbitals has 0 density, implying that the other two orbitals are involved in the edge states. BUT this does not include summing as in your PDOS equation... Is my analysis still correct? Jun 22 at 15:50
• Answering your last remark. Then you could do: w = wavefunction.conj() * wavefunction which would give you the PDOS weight for each orbital for that eigenstate. It seems that is all you want. Then you could do up = w[::3]; middle = w[1::3]; down = w[2::3] and then if np.allclose(middle, 0) you have a state which only lives on the up/down orbitals. So do you want to figure out which eigenvector is the edge state, or something different? Jun 22 at 16:13
• Regarding your first question if you both append evecs and eigs you'll get something like this: eigs = [nk, nbands] and evecs = [nk, number_of_orbitals, nbands], assuming it is the output directly from numpy/scipy. Then wavefunction = evecs[ik, :, iband]. Hope this makes sense? Jun 22 at 17:27
• Nice to see a lot of discussion. If the system prompts you to move to chat, I recommend this room. Please avoid clicking the button that creates a new room, since we'd really like to avoid this. Jun 22 at 18:50