Sorting energy bands in physics

Question

Obtaining the band structure is a standard problem in physics. Supposed there is a system which is described by a Hamiltonian matrix that depends on some system parameter $H(k)$. To find the energies of the system, one solves the Schrodinger's equation, which is just a eigenvalue equation: $$H(k)\Psi(k) = E(k)\Psi(k),$$ where $H$ is a Hermitian matrix, $\Psi(k)$ is the eigenfunction of $H$, and $E(k)$ is the eigenvalue (eigenenergy) of $H$.

If the Hamiltonian $H(k)$ is just a number, the plot $E(k)$ is a simple curve called dispersion. If $H$ is a $2\times 2$ matrix, $H(k)$ has two eigenvalues, $E_1(k)$ and $E_2(k)$. The plot $E$ vs $k$ containing the two curves is called the band structure of $H$.

It often happens that the bands $E_1(k)$ and $E_2(k)$ intersect for some values of k. After those intersections, it is not at all clear how to tell those bands apart. The situation gets even more complex for larger Hamiltonian matrices ($n\times n$).

There are some ways to tell the bands apart. One way is to check $\left<\Psi(k)\left| A \right|\Psi(k)\right>$ are continuous as a function of $k$, where A is the matrix of an operator defined in the same Hilbert space as $H$, which describes some property of the system. One can also check if the numerical derivatives of $E(k)$ on the left and right of the intersection points coincide.

But, these do not always give the correct result, and even very well established codes, the band sorting algorithms fail when i) there are too many bands, ii) when H is strongly coupled, i.e. it has large off-diagonal elements, comparable to the diagonal ones.

The question is, is there a better way to sort the bands, that will work in all situations, or making a finer $k$-grid is the only additional thing to be done?

Edit: how sorting works. Suppose we have the $k$-grid: $\{k_0, k_1, \dots, k_i, k_{i+1}\dots \}$. For each $k$, we diagonalize $H(k)$ and get two energies $E_1(k)$ and $E_2(k)$. Suppose $E_1(k)<E_2(k)$ for all $k_i$. Then no sorting is needed. First band will be the one consisting of all $E_1(k)$.

But, suppose that at $k_i$, $E_1(k_i) = E_2(k_i)+\Delta$, where $\Delta$ is just a small number. For $k_{i+1}$, just choosing the lower energy as the energy for the next band point is not sufficient. So, to sort bands, we imposed the condition that $E_1(k_i)$ is a continuous function and its derivatives are continuous. That works most of the time, but it is unclear to me what to do when it doesn't.

Answer 1

The open source Python package Kwant has built in routines to find the bands for any lattice Hamiltonian (tight-binding representation). It resolves the crossings by ensuring a continuous derivative. One developer of the package, Anton Akhmerov, has posting on this problem, which I found quite illuminating.

Answer 2

The algorithm used in Kwant, mentioned by Carlo uses adaptive sampling to resolve level crossings, and it it fairly involved. For for completeness I'll present a simpler approach that uses eigenvector continuity and works with the simplistic grid sampling.

The detailed description is in this blog post.

The key idea is to take a set of eigenvectors $|\psi_i(k)⟩$ at one $k$-point and figure out what is the most continuous way to connect it to the set of the eigenvectors at the neighboring $k$-point $|\psi_i(k + δk)⟩$. To do that we compute the matrix of overlaps $A_{ij} = |⟨\psi_i(k) | \psi_j(k + δk)⟩|^2$. We then find a bijection from $i$ to $j$ that maximizes the matched overlaps. This is an assignment problem, which is addressed by standard graph-theoretic algorithms.