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

Question

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

Answer 1

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. :)